#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <math.h>
#include <algorithm>
#include <time.h>

// GLS library 
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_monte_miser.h>


#define DEBUG
#define MODEL_upper
#define consumpEquiv_adj
#define MODEL_debt0
#define SOLVE_EGM
#define USE_stderror
#define MODEL_probs_educ5
#define MODEL_rational


#include "sandia_sparse.hpp"	        // For sparse grid

// Numerical Recipe 3rd: 
#include "nr3/nr3.h"
#include "nr3/gamma.h"
#include "nr3/gauss_wgts.h"


#ifdef USE_stderror

#include "nr3/ludcmp.h"
#include "nr3/savgol.h"
#include "nr3/fourier.h"
#include "nr3/convlv.h"

#endif


using namespace std;


#define Nvar 43 
#define Est_stepsize     0.02
#define SMM_replicates   30

#define NUM_measures_thetac  4   
#define NUM_measures_thetan  6   
#define NUM_unobs_loglike    2
#define NUM_icat_unobs       1   
#define NUM_measures_control 3  
#define Nvar_loglike            61  
#define Est_stepsize_loglike 	0.02 
#define Est_tol_loglike		    0.0001
#define Est_maxiter_loglike 	999000
#define NUM_loglike_ids         1605 
#define NUM_simsp_loglike       1000  




#define Est_tol		    0.0002
#define Est_maxiter 	6000
#define NUM_simsp		80  
				
		    	
// ... Guass-Hermite for price changes
#define neps_health      3   
		    	VecDoub xeps_health(neps_health);  
		    	VecDoub weps_health(neps_health);


// global variable used for Gauss-Laguerre integration integration over earnings shocks in solveModel_EGM.cpp; updated in getpar_XWEPS_logwf()    	
#define NUM_EPS_earnings  6   
		double XEPS_logwf[5][NUM_EPS_earnings]; 
		double WEPS_logwf[5][NUM_EPS_earnings]; 
		
// Gauss-Laguerre integration 
const int neps_earnings = NUM_EPS_earnings;  


// Grid for endogenous grid M 
#define NUM_MGrid  50  


#define NUM_DATA_init        1605
const int NUM_smmPerson  = SMM_replicates*NUM_DATA_init;    
#define NUM_smmcoef_wage     12
const int NUM_ismm = 107;  


#define Num_cf 			10 

//!!!!!!!! MODEL Primitives !!!!!!!!!
#define INIT_age        15 
#define NUM_age     	50 // 64 - 15 + 1 = 50 
#define NUM_age30       16 	
#define NUM_ageEduc     14 
const int NUM_age_ptr = NUM_ageEduc;      
#define NUM_ageMAXT     81 // 95 - INIT_age + 1


#define MIN_smmAge      0        

#define MAX_ageHS       20      
#define MIN_educ        7      
int MAX_initeduc  =  9;       


#define MAX_educ	    20			// maximum years of schooling
const int NUM_ieduc  = MAX_educ-MIN_educ+1;       
const double MAX_kstock = NUM_age;  

#define MAX_educpr      16
#define MIN_educpr      8   

#define MODEL_cgrid       500.0


const double MAX_savings[] = { 200000.0, 400000.0, 400000.0, 800000.0, 800000.0, };
#define MIN_savings      -100000.0      
#define SIZE_savingsStep        50000.0
#define SIZE_savingsStep_small  2000.0 
#define MAX_consumption   1.0e+30 
#define Factor_dUdC       1000.0


#define NUM_isavingspr  1
#define NUM_ieduccatpr  2 // 4


#define NUM_de         	2
#define DIM_theta       2       

#define NUM_dq          2       
#define NUM_educcat5    5

#define HOURS            1.0 // full-time hours
#define HOURSP           0.5 // part-time hours


//... Controls for Smolyak 
#define NUM_coefEmax  41  
#define DIM_smolyak   4  
#define RULE_smolyak  1    
int LEVEL_smolyak_vector[] = {2, 2, 2, 2};   

int ORDER_smolyak; int LEVEL_smolyak;  

double **GRID_smolyak; 

int **ChebyIndex; 
gsl_matrix * POLYINV_smolyak; 


const int NUM_networth_neg = 3;   
const int NUM_networth_pos = 5; 
const int NUM_networth = NUM_networth_pos + NUM_networth_neg + 1;    

double curv_networth_neg = 1.0; 
double curv_networth_pos = 2.0; 

#ifndef SOLVE_EGM
double coefEmax[NUM_age][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_networth][NUM_coefEmax]; 
#endif 

#ifdef MODEL_dEmax
double coefdEmax[NUM_age][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_networth][NUM_coefEmax]; 
#endif 

#ifdef MODEL_emaxQ0
double coefEmax_Q0[NUM_age][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_networth][NUM_coefEmax]; 
#endif 



double MIN_coefEmaxstate[NUM_age][NUM_ieduc][DIM_smolyak]; 
double MAX_coefEmaxstate[NUM_age][NUM_ieduc][DIM_smolyak];
double MAX_debt[NUM_age+1][NUM_ieduc], MAX_asset[NUM_age+1][NUM_ieduc];

double networthGrid[NUM_age+1][NUM_ieduc][NUM_networth];
double MAX_qstock[100];  

const double max_ci = 1.96;  


const double MIN_thetac = -2.596022;  
double       MAX_thetac =  -MIN_thetac;  
const double MIN_thetan = -2.001922;   
double       MAX_thetan = -MIN_thetan;  


double DATA_avecostdq[NUM_age+1][NUM_DATA_init]; double DATA_avetaxperpack[NUM_age+1][NUM_DATA_init]; 
int DATA_nlsyid[NUM_DATA_init]; 


int 	DATA_init_hgc[NUM_DATA_init], DATA_init_de_prev[NUM_DATA_init]; 
double  DATA_init_net_worth[NUM_DATA_init];
double  DATA_init_parents_hgc[NUM_DATA_init]; 
double	DATA_init_parents_net_worth[NUM_DATA_init]; 


double  DATA_init_measures_thetac [NUM_measures_thetac][NUM_DATA_init], 
 		DATA_init_measures_thetan [NUM_measures_thetan][NUM_DATA_init], 
        DATA_init_measures_control[NUM_measures_control][NUM_DATA_init], 
        DATA_init_icat_control    [NUM_icat_unobs][NUM_DATA_init]; 
double  DATA_init_icat_control_mean[NUM_icat_unobs]; 

double cutoff_measure_thetan[4]; 

double  DATA_famsize[100][NUM_educcat5]; 


const int NUM_age90 = 90-INIT_age+1; 
double DATA_deathrate[NUM_age90][5]; // CPS and Vital Statistics: 15 to 90: 90-15+1  
double DATA_deathrate_smoker_educ5[NUM_age90][5], DATA_deathrate_nonsmoker_educ5[NUM_age90][5]; 
double DATA_deathrate_smoker_educ3[NUM_age90][3], DATA_deathrate_nonsmoker_educ3[NUM_age90][3]; 
double DATA_deathrate_all[NUM_age90]; 
double Param_survival_new[13]; 



double TARGET_moments[NUM_ismm], TARGET_wgtmoments[NUM_ismm];   
	gsl_vector * olsY_logwage; 
	gsl_matrix * olsX_logwage;
	gsl_vector * olsCoef_logwage; 	

	gsl_vector * olsY_dqmore; 
	gsl_matrix * olsX_dqmore ; 
	gsl_vector * olsCoef_dqmore; 
	
	gsl_vector * olsY_savings	; 
	gsl_matrix * olsX_savings    ; 
	gsl_vector * olsCoef_savings; 
//	gsl_vector * olsY_savings_v2	; 
	gsl_matrix * olsX_savings_v2    ; 
	gsl_vector * olsCoef_savings_v2; 

int NUM_RetireT[] = {26, 26, 26, 26, 26}; 

//--------------------------------------------
//--- exogenous parameters
const double R_lend = 0.02; double R_borrow = 0.05;
const double R_interestRate = 0.05; 


double cbc_tc[]        = {5072.81859, 5072.81859, 10653.50853, 10653.50853, 11374, 11374, 11374, 11374};

double cbc_roomboard[] = {4539, 4539, 6532, 6532, 6532, 6532, 6532, 6532};  



const double cbc_grants =  1784.236*2.18056/1.666;  
double max_debtg_e[] = {2625.0, 3500.0, 5500.0, 5500.0, 10500.0, 10500.0, 10500.0, 10500.0};
double max_debtg[] = {23000.0, 138500.0-23000.0}; 
double maxpost_debtg = 0.0; 
double MAX_debt_gsl[NUM_ieduc]; 



const double c_chi = 7800; 


double cu_min; 

//------------------------------------------------------------------------------------
//--- structure parameters to be estimated: 
double delta0, delta1[3]; double delta2[2];
double cu_g; 
double cu_e_a; 
double cu_e[4], cu_e_rcost, cu_e_pr, cu_eh, cu_ke[2], cu_ka, cu_q[2], cu_qq, cu_q_alpha[2], cu_qstock;   
double cu_qe_cat; 
double cu_qh, cu_qde; 

double 	cu_de_dq, cu_e_qstock; 


double 	cu_e_alpha[DIM_theta]; 

double sigeps_ue, sigeps_ukp, sigeps_ukf, sigeps_uq;  

// ... parameter for the dynamics of log prices 
double dlogh    = 0.0417; 
double sigeps_h = 0.0877; 

const double MIN_logcprice = 5.9138-1.0; 
const double MAX_logcprice = 7.787573+1.0; 

double clogw_alpha[DIM_theta][4], clogw_skill[7]; 
double clogw_skill_graduate, clogw_ksquare_graduate; 

double clogw_beta[4]; double clogw_ndk; 

double sigeps_logwscale[5]; 
double sigeps_logwshape[] = {2.0, 2.0, 2.0, 2.0, 2.0};
const double MIN_epsw[]   = {1.0,  1.0,  1.0,  1.0,  1.0};


double cbc_ptr[24], sigeps_ptr; 
double cbc_gtr[4] ; 
double cdebt_beta[5], cdebt_alpha[2];  
double cdebt_min_age  =0.0; double cdebt_min_skill=0.0; double cdebt_min_0=0.0;



const int NUM_eps = 10; 
double  simp_rng_eps[NUM_simsp][NUM_eps]; 


double SMM_moments[NUM_ismm]; double SMM_sumf2; 

//!!!!!!!! Global Variables for Simulation of Individuals !!!!!!!!!!
double sm_rng_eps[NUM_age][NUM_smmPerson][NUM_eps]; 
double sm_rng_eps_unobs[NUM_smmPerson][NUM_unobs_loglike], sm_rng_eps_measure[NUM_age+1][NUM_smmPerson][NUM_unobs_loglike]; 

double sm_theta[DIM_theta][NUM_smmPerson];
double sm_wgt[NUM_smmPerson];
int sm_nlsy_num[NUM_smmPerson]; 


int sm_de_prev[NUM_age+1][NUM_smmPerson]; 			 
int sm_dq_prev[NUM_age+1][NUM_smmPerson]; 
		
int    sm_educ[NUM_age+1][NUM_smmPerson]; 
double sm_kstock[NUM_age+1][NUM_smmPerson], sm_expr[NUM_age+1][NUM_smmPerson];   
int sm_add[NUM_age+1][NUM_smmPerson]; 
double sm_qstock[NUM_age+1][NUM_smmPerson]; 		
double sm_cprice [NUM_age+1][NUM_smmPerson]; 
double sm_savings[NUM_age+1][NUM_smmPerson];

double sm_c[NUM_age][NUM_smmPerson];    


int	   sm_de[NUM_age][NUM_smmPerson]; 
double sm_dk[NUM_age][NUM_smmPerson];  
int    sm_dq[NUM_age][NUM_smmPerson]; 

double sm_wage[NUM_age][NUM_smmPerson];  
double sm_earnings[NUM_age][NUM_smmPerson];

double sm_probs[NUM_age][NUM_smmPerson]; 


int    sm_work_part[NUM_age][NUM_smmPerson], sm_work_full[NUM_age][NUM_smmPerson], sm_work_fullpart[NUM_age][NUM_smmPerson];  

double sm_skill[NUM_age][NUM_smmPerson], sm_debtlimit_p[NUM_age][NUM_smmPerson], sm_value[NUM_age][NUM_smmPerson]; 
double sm_epslogw[NUM_age][NUM_smmPerson], sm_hours[NUM_age][NUM_smmPerson]; 

int    sm_init_dq[NUM_smmPerson][2]; 


int    sm_isavingspr[NUM_smmPerson];
int    sm_ieduccatpr[NUM_smmPerson];
int    sm_educcat4_a30[NUM_smmPerson]; 

int    sm_educcat[2][4][NUM_smmPerson];  
int    sm_add0[NUM_age][NUM_smmPerson];

double sm_gtransfer[NUM_age][NUM_smmPerson], sm_gtrcmin[NUM_age][NUM_smmPerson], sm_ptransfer[NUM_age][NUM_smmPerson];

double sm_measure_thetac[NUM_smmPerson], sm_measure_thetan[NUM_smmPerson];
double sm_tax_cprice[NUM_age+1][NUM_smmPerson]; 
int sm_nlsyid[NUM_smmPerson]; 

double simp_rng_unobs_loglike[NUM_simsp_loglike][NUM_unobs_loglike]; 

double  cmu_measures_thetac[NUM_measures_thetac], 
		calpha_measures_thetac[NUM_measures_thetac], 
		csigma_error_thetac[NUM_measures_thetac]; 
double 	cmu_measures_thetan[NUM_measures_thetan], 
		calpha_measures_thetan[NUM_measures_thetan], 
		csigma_error_thetan[NUM_measures_thetan]; 
double  cmu_unobs[NUM_icat_unobs][NUM_unobs_loglike]; 
double  cmu_unobs0[NUM_unobs_loglike]; 

const int ncov = (NUM_unobs_loglike*NUM_unobs_loglike-NUM_unobs_loglike)/2; 
double  clvar_unobs[NUM_unobs_loglike]; 
double  clcov_unobs[ncov]; 
gsl_matrix * MAT_lcov_unobs     = gsl_matrix_alloc(NUM_unobs_loglike, NUM_unobs_loglike); 
gsl_matrix * MAT_cov_unobs      = gsl_matrix_alloc(NUM_unobs_loglike, NUM_unobs_loglike); 


double cbeta_measures_thetac[NUM_measures_thetac][NUM_measures_control]; 
double cbeta_measures_thetan[NUM_measures_thetan][NUM_measures_control];  



	const int NUM_res = NUM_DATA_init*NUM_simsp_loglike;	
	double measure_hat_c[NUM_measures_thetac][NUM_res]; 
	double measure_hat_n[NUM_measures_thetan][NUM_res]; 
	double measure_hat2_n[NUM_measures_thetan][NUM_res]; 
	



#ifdef USE_mpi
int mpi_nproc, mpi_nii, mpi_itop, mpi_iend;
#else 
const int mpi_id = 0;
#endif

double cc_taxPolicy_incomeTest = 1.0e+30;  
int cc_taxPolicy = 0; 
double gtr_sub = 0.0; double gtr_csub = 0.0; double cbc_sub = 0.0;  
double TAXrate_y = 0.0; double TAXrate_dq = 0.0; 


double 	everAVEdq_a15[Num_cf], everAVEdq_a20[Num_cf], everAVEdq_a25[Num_cf], everAVEdq_a30[Num_cf]; 


	double AVEagedeath_a15[Num_cf], AVEyearsleft_a15[Num_cf], everAVEagedeath_a15[Num_cf], everAVEyearsleft_a15[Num_cf]; 

	double AVEdqever_a25[Num_cf]; 
	double AVEdq_a15[Num_cf]; 
	double AVEdq_a20[Num_cf]; 
	double AVEdq_a25[Num_cf]; 
	double AVEdq_a30[Num_cf]; 

	double AVEhgc_a15[Num_cf]; 
	double AVEhgc_a20[Num_cf]; 
	double AVEhgc_a25[Num_cf]; 
	double AVEhgc_a30[Num_cf]; 
	

	double AVEc_a15[Num_cf]; 
	double AVEc_a20[Num_cf]; 
	double AVEc_a25[Num_cf]; 
	double AVEc_a30[Num_cf]; 


	double AVEnetworth_a15[Num_cf]; 
	double AVEnetworth_a20[Num_cf]; 
	double AVEnetworth_a25[Num_cf]; 
	double AVEnetworth_a30[Num_cf]; 

	double AVEearnings_a15[Num_cf]; 
	double AVEearnings_a20[Num_cf]; 
	double AVEearnings_a25[Num_cf]; 
	double AVEearnings_a30[Num_cf]; 


	double AVEvalue_a15[Num_cf]; 
	double AVEvalue_a20[Num_cf]; 
	double AVEvalue_a25[Num_cf]; 
	double AVEvalue_a30[Num_cf]; 


	double AVEqstock_a15[Num_cf]; 
	double AVEqstock_a20[Num_cf]; 
	double AVEqstock_a25[Num_cf]; 
	double AVEqstock_a30[Num_cf]; 
	
	
	double AVEadd_a15[Num_cf]; 
	double AVEadd_a20[Num_cf]; 
	double AVEadd_a25[Num_cf]; 
	double AVEadd_a30[Num_cf]; 
	
double VARlogearnings[Num_cf], VARlogskill[Num_cf], VARlogc[Num_cf], VARloghgc[Num_cf], Probhgc12_a19[Num_cf], Probhgc13_a21[Num_cf], Probhgc16_a25[Num_cf];
double AVEearnings[Num_cf], AVEc[Num_cf], AVEskill[Num_cf], AVEhgc[Num_cf], AVEnetworth[Num_cf]; 
double GiniWealth[Num_cf], GiniEarnings[Num_cf]; 
double AVEinitvalue[Num_cf], VARloginitvalue[Num_cf];  double AVEvaluea30[Num_cf], VARlogvaluea30[Num_cf]; 
double WorkpartEnroll[Num_cf]; 

double VARlogcprice[Num_cf], VARlogqstock[Num_cf];
double AVEhealth[Num_cf], AVEqstock[Num_cf], AVEdq[Num_cf], AVEinitdq[Num_cf], AVEeverdq[Num_cf], AVEdq1630[Num_cf]; 

double  AVEhgc11less[Num_cf], AVEhgc12[Num_cf], AVEhgc1315[Num_cf], AVEhgc16more[Num_cf], AVEhgc11less_everdq[Num_cf], AVEhgc12_everdq[Num_cf], AVEhgc1315_everdq[Num_cf], AVEhgc16more_everdq[Num_cf], AVEhgc12more[Num_cf], AVEhgc12more_everdq[Num_cf]; 
      

double everProbhgc12_a19[Num_cf], everProbhgc13_a21[Num_cf], everProbhgc16_a25[Num_cf], everAVEeverdq[Num_cf], everAVEdq1630[Num_cf], everAVEdq[Num_cf], everAVEearnings[Num_cf], everAVEc[Num_cf], everAVEnetworth[Num_cf], everAVEhgc[Num_cf], everAVEadd[Num_cf], everAVEqstock[Num_cf], everAVEinitvalue[Num_cf], 
everAVEvaluea30[Num_cf], everAVEskill[Num_cf]; 

double dqa30Probhgc12_a19[Num_cf], dqa30Probhgc13_a21[Num_cf], dqa30Probhgc16_a25[Num_cf], dqa30AVEdqa30dq[Num_cf], dqa30AVEdq1630[Num_cf], dqa30AVEdq[Num_cf], dqa30AVEearnings[Num_cf], dqa30AVEc[Num_cf], dqa30AVEnetworth[Num_cf], dqa30AVEhgc[Num_cf], dqa30AVEadd[Num_cf], dqa30AVEqstock[Num_cf], dqa30AVEinitvalue[Num_cf], 
dqa30AVEvaluea30[Num_cf], dqa30AVEskill[Num_cf]; 

double typecatAVEhgc[Num_cf][2], typecatAVEskill[Num_cf][2], typecatAVEqstock[Num_cf][2], typecatAVEdq[Num_cf][2], typecatAVEearnings[Num_cf][2], typecatAVEc[Num_cf][2], typecatAVEnetworth[Num_cf][2], typecatAVEinitvalue[Num_cf][2], typecatAVEvaluea30[Num_cf][2]; 
double typecatProbhgc12_a19[Num_cf][2], typecatProbhgc13_a21[Num_cf][2], typecatProbhgc16_a25[Num_cf][2]; 

double educcatAVEhealth[Num_cf][4], educcatAVEskill[Num_cf][4], educcatAVEqstock[Num_cf][4], educcatAVEdq[Num_cf][4], educcatAVEearnings[Num_cf][4], educcatAVEc[Num_cf][4], educcatAVEnetworth[Num_cf][4], educcatAVEinitvalue[Num_cf][4], educcatAVEvaluea30[Num_cf][4]; 
double AVEeduccat4[Num_cf][4];
double educcatAVElogcprice[Num_cf][4], educcatAVElogearnings[Num_cf][4], educcatAVElogc[Num_cf][4];


double logcprice_investment(int educ, int dq, int de, double earnings, int ia); 
double get_measure_healthstat(double logcprice, double rng, double measures_control[]); 

double get_kstock(int educ, int ia); 
double get_skill(int educ, double theta[2], int ia, double epsw0); 
double get_debtlimit_endogenous_maxcoef(const int educ_next1, double theta[2], const int ia); 
double get_debtlimit_ptotal(double theta[2], const int ia,  const int educ, const int de, const int dq, const double qstock); 


void readpar(const int n, ifstream &infile);
void printpar(const char *fname); 
double get_appEmax_savingsR(int de_prev0, int isavingspr0, int ieduccatpr0, int ia, int ieduc, double theta[], double kstock, double savingsR, double logcprice, double qstock); 
double get_appEmax(int de_prev0, int isavingspr0, int ieduccatpr0, int ia, int ieduc, double theta[], double kstock, double savings, double logcprice, double qstock); 
double get_dEmax_dS(int ia, int ieduc, double theta[], double kstock, double savings, double *polycoef); 


double get_stateSMOLYAK(const int iss, double theta[], const double logcprice, const double qstock, const int ia, const int ieduc); 
int get_stateMODEL(const double state[], double &thetac, double &thetan, double &logcprice, double &qstock, const int ia, const int ieduc);
void get_coefEmaxstate();  

void get_coefEmax(); 
void get_smOpt(int miniper, int maxiper, double coefemaxv[NUM_age][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_networth][NUM_coefEmax]); 
double smm_obj(const gsl_vector *x, void *void_params);
double get_sumf2(); 
int taxPolicy_model (); 
double get_dqtax_rate(int ia, int dq); 
double get_incomeTax(double inc_wages, double inc_s); 
void write_coefEmax(const char *fname); 
inline double get_savings_BY_savings_norm(const double savings_norm, const int dummy_negSavings, int ia, int ieduc); 
inline double get_rng_epsw(int educ, double rng_eps);
double get_prob_myopic(double theta[], int hgc0, int de_prev0, int ieduccatpr0); 
double get_lifeExpectancy_static(const int ia, const int educ, const double qstock, const int dq, const int de); 
double get_survivalProb(int ia, int educ0, double qstock, int dq, int de); 

#include "myFunction.cpp"	       	   // function for standard errors
#include "mySmolyak.cpp"

#include "loglike_measures.cpp"	       // measurement system



struct struct_infoset{
// health is the price variable. 	
	
    double *theta; 
    double health, logcprice;  
	
	double delta;     
    
    
	int ia, de_prev; 
	
    int educ; 

    double kstock;

    double qstock;   
	double skill; 	
	
	double savings; 


	int isavingspr; 	    
	int ieduccatpr; 
	
	
	double eps_ue, eps_uq;    
	
	double eps_logwf;
	
	double eps_health_next; 
	
	double eps_ptr, eps_ptrdummy;  
	
	double debtlimitLag; 
	
}; 


struct struct_foc{
	
		
	double cashathand; 
			
	int de, dq; double dk;  
	
	struct_infoset infoset; 
	
	int isavingspr_next, ieduccatpr_next; 
	
}; 



void printpar(const char *fname){

	FILE * outfile; outfile  = fopen (fname, "w"); 
        
    fprintf(outfile, "%.8f \n", cu_min);

	fprintf(outfile, "%.8f \n", cu_g);
		
	fprintf(outfile, "%.8f \n", clogw_skill_graduate); 
	fprintf(outfile, "%.8f \n", clogw_ksquare_graduate); 
	fprintf(outfile, "%.8f \n", delta2[0]);
	fprintf(outfile, "%.8f \n", delta2[1]);

	
	fprintf(outfile, "%.8f \n", delta0);

	for (int idelta=0; idelta<2; idelta++){
		fprintf(outfile, "%.8f \n", delta1[idelta]);	
	}
	

			
	fprintf(outfile, "%.8f \n", cu_e[0]);	
	fprintf(outfile, "%.8f \n", cu_e[1]);
	fprintf(outfile, "%.8f \n", cu_e[2]);
    fprintf(outfile, "%.8f \n", cu_e[3]);
	fprintf(outfile, "%.8f \n", cu_e_rcost);
	fprintf(outfile, "%.8f \n", cu_e_alpha[0]);	
	fprintf(outfile, "%.8f \n", cu_e_alpha[1]);
	
	fprintf(outfile, "%.8f \n", cu_e_pr);	

	fprintf(outfile, "%.8f \n", cu_e_a);	
	

	fprintf(outfile, "%.8f \n", cu_ka);

	
	fprintf(outfile, "%.8f \n", cu_ke[0]);	
	fprintf(outfile, "%.8f \n", cu_ke[1]);	
			
	fprintf(outfile, "%.8f \n", cu_q[0]);	
	fprintf(outfile, "%.8f \n", cu_q[1]);	
		
	fprintf(outfile, "%.8f \n", cu_qe_cat);
			
	fprintf(outfile, "%.8f \n", cu_qq);		
	fprintf(outfile, "%.8f \n", cu_q_alpha[0]);	
	fprintf(outfile, "%.8f \n", cu_q_alpha[1]);	
	
	
	fprintf(outfile, "%.8f \n", sigeps_ue);	
	fprintf(outfile, "%.8f \n", sigeps_uq);	

	//... wage equation
	for (int iw=0; iw<7; iw++){
	fprintf(outfile, "%.8f \n", clogw_skill[iw] );	
	}

	//... wage equation
	fprintf(outfile, "%.8f \n", clogw_beta[0]);
	fprintf(outfile, "%.8f \n", clogw_beta[1]);	
	fprintf(outfile, "%.8f \n", clogw_beta[2]);	
	fprintf(outfile, "%.8f \n", clogw_beta[3]);	
	
	
	fprintf(outfile, "%.8f \n", clogw_ndk);	
	
	
	for (int ieduccat = 0; ieduccat < 4; ieduccat ++) {
			
	fprintf(outfile, "%.8f \n", clogw_alpha[0][ieduccat]);	
	fprintf(outfile, "%.8f \n", clogw_alpha[1][ieduccat]);
		
	}
	
	for (int ieduccat = 0; ieduccat < 5; ieduccat ++) {
			
	fprintf(outfile, "%.8f \n", sigeps_logwscale[ieduccat]);	
	
	}	

	for (int ieduccat = 0; ieduccat < 5; ieduccat ++) {
			
	fprintf(outfile, "%.8f \n", sigeps_logwshape[ieduccat]);	
	
	}
    
	
    fprintf(outfile, "%.8f \n", cu_qstock);


	fprintf(outfile, "%.8f \n", cu_de_dq);
	fprintf(outfile, "%.8f \n", cu_e_qstock);	 

	 
	fclose(outfile); 
	
}


void readpar(const int n, ifstream &infile){

	double temp_x[n]; 
	
    for(int i=0; i< n && !infile.eof(); i++) infile >> temp_x[i];
    
	int i = 0;
            
    cu_min = temp_x[i];         i++;   

	cu_g    = temp_x[i]; 		i++; 
	
	clogw_skill_graduate = 	temp_x[i]; 		i++;	
	clogw_ksquare_graduate = temp_x[i]; 		i++;			
	  
	delta2[0] = temp_x[i]; 		i++;
	delta2[1] = temp_x[i]; 		i++;
	
	delta0 = temp_x[i]; 		i++; 
	
	for (int idelta=0; idelta<2; idelta++){
		delta1[idelta] = temp_x[i];	i++; 
	}
	
		
	cu_e[0] = temp_x[i]; 		i++;       	
	cu_e[1] = temp_x[i]; 		i++;       	
	cu_e[2] = temp_x[i]; 		i++;       	
    cu_e[3] = temp_x[i];        i++;
    
	cu_e_rcost = temp_x[i]; 		i++; 		
	cu_e_alpha[0] = temp_x[i]; 		i++; 	
	cu_e_alpha[1] = temp_x[i]; 		i++; 	
	
	
	cu_e_pr   = temp_x[i]; 		i++; 	
	
	cu_e_a  = temp_x[i]; 		i++; 
	

	cu_ka = temp_x[i]; 		i++; 		
	

	cu_ke[0] = temp_x[i]; 		i++; 
	cu_ke[1] = temp_x[i]; 		i++; 
	

	
	cu_q[0] = temp_x[i]; 		i++; 		
	cu_q[1] = temp_x[i]; 		i++; 		
	
	cu_qe_cat  = temp_x[i]; 		i++; 		
	
	cu_qq   = temp_x[i]; 		i++; 		
	cu_q_alpha[0] = temp_x[i]; 		i++; 
	cu_q_alpha[1] = temp_x[i]; 		i++; 
	
	
	sigeps_ue = temp_x[i]; 		i++; 
	sigeps_uq = temp_x[i]; 		i++; 

	
	for (int iw=0; iw<7; iw++){
	clogw_skill[iw]  = temp_x[i]; 		i++;   
	} 

	
	clogw_beta[0]    = temp_x[i]; 		i++;    
	clogw_beta[1]    = temp_x[i]; 		i++;   
	clogw_beta[2]    = temp_x[i]; 				i++;    
	clogw_beta[3]    = temp_x[i]; 				i++;  
	
	
	clogw_ndk  = temp_x[i]; 		i++; 
	
	
	for (int ieduccat = 0; ieduccat < 4; ieduccat ++) {
	
	clogw_alpha[0][ieduccat] = temp_x[i]; 		i++; 	
	clogw_alpha[1][ieduccat] = temp_x[i]; 		i++; 	
		
	}

	for (int ieduccat = 0; ieduccat < 5; ieduccat ++) {
	
	sigeps_logwscale[ieduccat]    = temp_x[i]; 		i++;  

	}
	
	for (int ieduccat = 0; ieduccat < 5; ieduccat ++) {
	
	sigeps_logwshape[ieduccat]    = temp_x[i]; 		i++;  

	}
	
    cu_qstock = temp_x[i];         i++;

	cu_de_dq    = temp_x[i];         i++;
	cu_e_qstock = temp_x[i];         i++;


	
	#ifdef DEBUG
	
	printf("\n\n readpar: n = %3d \n\n", i); 
	
	#endif
	
}





double get_logcprice_next(int ia, int educ, int dq, int de, double dk, double skill, double theta[], double logcprice, double rng){
    
  double logcprice_next = logcprice + dlogh + rng * sigeps_h; 
  
  if ( logcprice_next < MIN_logcprice ) {
		logcprice_next = MIN_logcprice; 
  } else if (logcprice_next > MAX_logcprice){
  		logcprice_next = MAX_logcprice;
  }
  
  return logcprice_next; 

}

inline double get_eps_health_next(double logcprice, double logcprice_next){

  return (logcprice_next - (logcprice + dlogh) )/ sigeps_h; 

}

double get_logcprice_next(double logcprice, double rng){
    
  double logcprice_next = logcprice + dlogh + rng * sigeps_h; 
  
  if ( logcprice_next < MIN_logcprice ) {
		logcprice_next = MIN_logcprice; 
  } else if (logcprice_next > MAX_logcprice){
  		logcprice_next = MAX_logcprice;
  }
  
  return logcprice_next; 

}



void init_par(){

	cc_taxPolicy = 0; 	
	TAXrate_y = 0.0; 
	TAXrate_dq = 0.0;
	gtr_sub = 0.0; gtr_csub = 0.0; cbc_sub =0.0; 	
	
	cout << "\n\n Parameter value initialization completed \n\n"; 
	
}


void initpar_loglike(){

	gsl_rng * rng; 	gsl_rng_env_setup();
	rng = gsl_rng_alloc (gsl_rng_taus);
	int seed = 123;		gsl_rng_set( rng, seed ) ; 
	for (int issp=0; issp < NUM_simsp_loglike; issp++){
 		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
 	 		simp_rng_unobs_loglike[issp][iunobs] = gsl_ran_ugaussian(rng); 	 	
 	 	}
 	 }
	
		
	for(int im = 0; im < NUM_measures_thetac; im++){
		
		cmu_measures_thetac[im] = 0.0; 
		
		// factor loading in the measurement equations
		calpha_measures_thetac[im] = 1.0; 
		
		// measurement errors
		csigma_error_thetac[im] = 1.0; 
		
		// covariance of thetac 
		for (int ic = 0; ic < NUM_measures_control; ic++){
			cbeta_measures_thetac[im][ic] = 0.0; 
		}
	}
	
	for(int im = 0; im < NUM_measures_thetan; im++){
		
		cmu_measures_thetan[im] = 0.0; 
		
		// factor loading in the measurement equations
		calpha_measures_thetan[im] = -1.0; 

		// measurement errors
		csigma_error_thetan[im] = 1.0; 
		
		// covariance of thetan 
		for (int ic = 0; ic < NUM_measures_control; ic++){
			cbeta_measures_thetan[im][ic] = 0.0; 
		}
	}
	
	

	
	// unconditional distribution of unobservables
	for (int icat =0; icat < NUM_icat_unobs; icat++){
		
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
			cmu_unobs[icat][iunobs] = 0.0; 
		}
	}

	gsl_matrix_set_zero (MAT_cov_unobs);
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++)	gsl_matrix_set(MAT_cov_unobs, iunobs, iunobs, 1.0);


	gsl_matrix_memcpy (MAT_lcov_unobs, MAT_cov_unobs); 
	gsl_linalg_cholesky_decomp (MAT_lcov_unobs); 
	for (int iunobs=0; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol = NUM_unobs_loglike-1; icol > iunobs ; icol--){
			gsl_matrix_set( MAT_lcov_unobs, iunobs, icol, 0.0); 	// ... set upper triangular matrix to zero 
		}
	}
	
	int icount =0; 
	
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
		clvar_unobs[iunobs] = gsl_matrix_get( MAT_lcov_unobs, iunobs, iunobs); 
		
	}
	icount =0; 
	for (int iunobs=1; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol =0; icol < iunobs ; icol++){
			clcov_unobs[icount] = gsl_matrix_get( MAT_lcov_unobs, iunobs, icol); 
			icount++; 
		}
	}
	

}


void init_data(){

    int n; int icnt; 
    ifstream infile;   
    
  
    n = NUM_DATA_init; 	
	
	
	// The init_distribution_NLSY97.txt contains BLS restricted access data and is not available for distribution 
	infile.open("out_data_new/init_distribution_NLSY97.txt"); 	
    if(!infile) {cerr<<"\n Restricted access data file (init_distribution_NLSY97.txt) is not available! \n\n"; exit(1); }

	for (int iobs=0; iobs < n * NUM_age && !infile.eof(); iobs++){	
		
		int cnt, nlsyid, age; 
		
		infile >> cnt >> nlsyid >> age ; 
		
		const int i = cnt - 1; 
		
		DATA_nlsyid[i] = nlsyid; 
		
		
		int it = age - INIT_age; 		
		infile >> DATA_avecostdq[it][i] >> DATA_avetaxperpack[it][i]; 

		
		double measure_thetac[NUM_measures_thetac]; 
		for (int ivec=0; ivec < NUM_measures_thetac; ivec++){
			
			infile >> measure_thetac[ivec]; 
			
			if (age == INIT_age) { DATA_init_measures_thetac[ivec][i] = measure_thetac[ivec]; 
	
			} 
		}
		

		
		double measure_thetan[NUM_measures_thetan]; 
		for (int ivec=0; ivec < NUM_measures_thetan; ivec++){

			infile >> measure_thetan[ivec]; 

			if (age == INIT_age) { DATA_init_measures_thetan[ivec][i] = measure_thetan[ivec]; 
			}

		}


				
		double measures_control[NUM_measures_control]; 
		for (int ivec=0; ivec < NUM_measures_control; ivec++){
			
			infile >> measures_control[ivec]; 
			
			if (age == INIT_age) { DATA_init_measures_control[ivec][i] = measures_control[ivec]; 
			}	
		}
		

			
		double icatcontrol[NUM_icat_unobs]; 	
		for (int ivec=0; ivec < NUM_icat_unobs; ivec++){
			infile >> icatcontrol[ivec]; 
			if (age==INIT_age) { DATA_init_icat_control[ivec][i] = icatcontrol[ivec]; 			
			}
		}
		

					
		int hgc; double parents_hgc; 
		infile >>  parents_hgc >>  hgc ;
		if (age==INIT_age) {
			DATA_init_parents_hgc[i]       = parents_hgc;		
			DATA_init_hgc[i]               = hgc;
		}

	}	    
	//... fill the future prices for the purpose of simulation at later part of lifecycle 
	for (int i=0; i < n; i++){	
 		for (int it=33-INIT_age+1; it<NUM_age+1; it++){
 		DATA_avecostdq[it][i]     = DATA_avecostdq[33-INIT_age][i]; 
 		DATA_avetaxperpack[it][i] = DATA_avetaxperpack[33-INIT_age][i]; 
 		}
		
		DATA_init_parents_net_worth[i] = 2; 
		DATA_init_net_worth[i] = 0.0; 
		DATA_init_de_prev[i] =  1; 
	}
	infile.close(); 
	
	
	for (int ivec=0; ivec < NUM_measures_control; ivec++){
		double m = my_mean(NUM_DATA_init, DATA_init_measures_control[ivec]); 
		if (fabs(m) > 0.02) { printf("\n\n Warning: control variables in measurement equation are not de-meaned: avg(%d) = %f ... \n exiting ... \n\n", ivec, m); exit(1); }  
	}	
		
	for (int icat=0; icat<NUM_icat_unobs; icat++){
		DATA_init_icat_control_mean[icat] = my_mean(NUM_DATA_init, DATA_init_icat_control[icat]); 
	}    

	
	const int neduc_available = 5; 
		
	double DATA_totaldeath[NUM_age90][neduc_available], DATA_totalpop[NUM_age90][neduc_available], DATA_totaldeath_nonsmoker[NUM_age90][neduc_available]; 
	for(int ia=0; ia < NUM_age90; ia++){
	DATA_deathrate_all[ia] = 0.0; 
	for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
		DATA_totaldeath[ia][ieduccat] = 0.0; 
		DATA_totalpop[ia][ieduccat] = 0.0; 
		DATA_totaldeath_nonsmoker[ia][ieduccat] = 0.0; 
	}}
	
	
	// 2000 census age 90 population is based on (age 89 population from census 2000) - (age 90 death from mort2001): (it does not matter because HRS is used for 50+ )	
	DATA_totalpop[NUM_age90-1][0] = 48786 - 7097.0; 
	DATA_totalpop[NUM_age90-1][1] = 40666 - 5505.0; 
	DATA_totalpop[NUM_age90-1][2] = 11463 - 1721.0; 
	DATA_totalpop[NUM_age90-1][3] = 8731 - 1483.0; 
	DATA_totalpop[NUM_age90-1][4] = 8428 - 1187.0; 	
	
	// start from age 15 
	infile.open("out_data_new/usa2000_Totals_male_educcat5_age15-89.txt", ifstream::in); 
    if(!infile) {cerr<<"\n Error reading the file usa2000_Totals_male_educcat3_age15-89 \n"; exit(1); }
	for(int i=0; i < NUM_age90*neduc_available && !infile.eof(); i++){
		
		int age; int educcat3; double pop; 
		infile >> age >> educcat3 >> pop; 

		int ieduccat = educcat3 -1; 
		int ia = age-INIT_age; 		
		DATA_totalpop[ia][ieduccat] = pop; 
		
	//	printf("\n %d, age %3d, ieduccat %d, pop %.2f ", i, ia+INIT_age, ieduccat, DATA_totalpop[ia][ieduccat]); 

	}
	infile.close();
		
	
	// start from age 16 
	infile.open("out_data_new/mort2001_Totals_male_educcat5_age16-91.txt", ifstream::in); 
    if(!infile) {cerr<<"\n Error reading the file mort2001_Totals_male_educcat5.txt \n"; exit(1); }
	for(int i=0; i <  NUM_age90*neduc_available && !infile.eof(); i++){
		
		int age; int educcat5; double alldeath; double cancerdeath; 
		infile >> age >> educcat5 >> alldeath >> cancerdeath; 
		
		int ieduccat = educcat5 -1; 
		int ia = age-(INIT_age+1); 
		
		DATA_totaldeath[ia][ieduccat] = alldeath; // 
		
		DATA_totaldeath_nonsmoker[ia][ieduccat] = alldeath - cancerdeath; 
		
	}

	double DATA_totaldeath_educ3[NUM_age90][3], DATA_totalpop_educ3[NUM_age90][3], DATA_totaldeath_nonsmoker_educ3[NUM_age90][3]; 
	for(int ia=0; ia < NUM_age90-1; ia++){
				
				DATA_totaldeath_nonsmoker_educ3[ia][0] =  DATA_totaldeath_nonsmoker[ia][0]; 
				DATA_totaldeath_educ3[ia][0]   = DATA_totaldeath[ia][0]; 
				DATA_totalpop_educ3[ia][0]     = DATA_totalpop[ia][0]; 
				
				DATA_totaldeath_nonsmoker_educ3[ia][1] =  DATA_totaldeath_nonsmoker[ia][1] + DATA_totaldeath_nonsmoker[ia][2]; 
				DATA_totaldeath_educ3[ia][1]   = DATA_totaldeath[ia][1] + DATA_totaldeath[ia][2]; 
				DATA_totalpop_educ3[ia][1]     = DATA_totalpop[ia][1] + DATA_totalpop[ia][2]; 

				DATA_totaldeath_nonsmoker_educ3[ia][2] =  DATA_totaldeath_nonsmoker[ia][3] + DATA_totaldeath_nonsmoker[ia][4]; 
				DATA_totaldeath_educ3[ia][2]   = DATA_totaldeath[ia][3] + DATA_totaldeath[ia][4]; 
				DATA_totalpop_educ3[ia][2]     = DATA_totalpop[ia][3] + DATA_totalpop[ia][4]; 
				
	}
		
	// ... assign age 90 death rate to be 1. 
	for (int ieduccat =0; ieduccat<5; ieduccat++) {
		DATA_deathrate[NUM_age90-1][ieduccat] = 1.0; 
		DATA_deathrate_smoker_educ5[NUM_age90-1][ieduccat] = 1.0; 
		DATA_deathrate_nonsmoker_educ5[NUM_age90-1][ieduccat] = 1.0; 
	}
	for(int ia=0; ia < NUM_age90-1; ia++){
		
		
		double totaldeath = 0.0; double totalpop = 0.0; 
		int current_neduc; 
		
		
		current_neduc = 1 + (MAX_initeduc+ia >= 12) + (MAX_initeduc+ia >= 13) + (MAX_initeduc+ia >= 16) + (MAX_initeduc+ia >= 17);
		totaldeath = 0.0; totalpop = 0.0; 
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			totaldeath += DATA_totaldeath[ia][ieduccat]; 
			totalpop   += DATA_totalpop[ia][ieduccat]; 

			if ( ieduccat <  current_neduc ){
				
		
				if (DATA_totaldeath[ia][ieduccat] > DATA_totalpop[ia][ieduccat] || DATA_totalpop[ia][ieduccat] < 1000.0){printf("\n Error !!! check death data: ia %d, ieduccat %d, death %f, pop %f !!!\n exiting \n", ia, ieduccat,DATA_totaldeath[ia][ieduccat], DATA_totalpop[ia][ieduccat] ); exit(1); }
				
				DATA_deathrate_smoker_educ5[ia][ieduccat]    = DATA_totaldeath[ia][ieduccat]/DATA_totalpop[ia][ieduccat]; 
				DATA_deathrate_nonsmoker_educ5[ia][ieduccat] = DATA_totaldeath_nonsmoker[ia][ieduccat]/DATA_totalpop[ia][ieduccat]; 
				
			} else {
				
				// assign to the nearest deathrate if education category is not available 
				
				if (ieduccat == 0){printf("\n Errors in setup.cc \n exiting \n"); exit(1); }
				
				DATA_deathrate_smoker_educ5[ia][ieduccat]    = DATA_deathrate_smoker_educ5[ia][ieduccat-1]; 
				DATA_deathrate_nonsmoker_educ5[ia][ieduccat] = DATA_deathrate_nonsmoker_educ5[ia][ieduccat-1]; 
				
			}
			
			if (DATA_deathrate_smoker_educ5[ia][ieduccat] < DATA_deathrate_nonsmoker_educ5[ia][ieduccat] ){printf("\n\n Error!! Please check the data: DATA_deathrate_smoker_educ5[%d][%d] %f < DATA_deathrate_nonsmoker_educ5[ia][ieduccat] %f \n\n ", ia, ieduccat, DATA_deathrate_smoker_educ5[ia][ieduccat], DATA_deathrate_nonsmoker_educ5[ia][ieduccat] ); exit(1); }

		} // ieduccat
		DATA_deathrate_all[ia] = totaldeath/totalpop; 		
		
		current_neduc = 1 + (MAX_initeduc+ia >= 12) + (MAX_initeduc+ia >= 16) ;
		totaldeath = 0.0; totalpop = 0.0; 
		for (int ieduc3 =0; ieduc3<3; ieduc3++){
			totaldeath += DATA_totaldeath[ia][ieduc3]; 
			totalpop   += DATA_totalpop[ia][ieduc3]; 
			
			if ( ieduc3 <  current_neduc ){
				
				if (DATA_totaldeath_educ3[ia][ieduc3] > DATA_totalpop_educ3[ia][ieduc3] || DATA_totalpop_educ3[ia][ieduc3] < 1000.0){printf("\n Error !!! check death data: ia %d, ieduc3 %d, death %f, pop %f !!!\n exiting \n", ia, ieduc3, DATA_totaldeath_educ3[ia][ieduc3], DATA_totalpop_educ3[ia][ieduc3] ); exit(1); }
				
				DATA_deathrate_smoker_educ3[ia][ieduc3]    = DATA_totaldeath_educ3[ia][ieduc3]/DATA_totalpop_educ3[ia][ieduc3]; 
				DATA_deathrate_nonsmoker_educ3[ia][ieduc3] = DATA_totaldeath_nonsmoker_educ3[ia][ieduc3]/DATA_totalpop_educ3[ia][ieduc3]; 
				
			} else {
				
				// assign to the nearest deathrate if education category is not available 
				
				if (ieduc3 == 0){printf("\n Errors in setup.cc \n exiting \n"); exit(1); }
				
				DATA_deathrate_smoker_educ3[ia][ieduc3]    = DATA_deathrate_smoker_educ3[ia][ieduc3-1]; 
				DATA_deathrate_nonsmoker_educ3[ia][ieduc3] = DATA_deathrate_nonsmoker_educ3[ia][ieduc3-1]; 
				
			}

		
		}
		

		if (ia+INIT_age<= 50){
		
		  for (int ieduc3 =0; ieduc3<2; ieduc3++){
		
			if (DATA_deathrate_smoker_educ3[ia][ieduc3+1] > DATA_deathrate_smoker_educ3[ia][ieduc3] ){
				DATA_deathrate_smoker_educ3[ia][ieduc3+1] = DATA_deathrate_smoker_educ3[ia][ieduc3]; 
			}

			if (DATA_deathrate_nonsmoker_educ3[ia][ieduc3+1] > DATA_deathrate_nonsmoker_educ3[ia][ieduc3] ){
				DATA_deathrate_nonsmoker_educ3[ia][ieduc3+1] = DATA_deathrate_nonsmoker_educ3[ia][ieduc3]; 
			}
					
		  }
		  
		  
		  
		  
		} // ia+INIT_age<= 50
		
		
		if (ia+INIT_age<= 50){
		
		  for (int ieduccat =0; ieduccat<5; ieduccat++){
		
			if (DATA_deathrate_smoker_educ5[ia][ieduccat+1] > DATA_deathrate_smoker_educ5[ia][ieduccat] ){
				DATA_deathrate_smoker_educ5[ia][ieduccat+1] = DATA_deathrate_smoker_educ5[ia][ieduccat]; 
			}
			
			if (DATA_deathrate_nonsmoker_educ5[ia][ieduccat+1] > DATA_deathrate_nonsmoker_educ5[ia][ieduccat] ){
				DATA_deathrate_nonsmoker_educ5[ia][ieduccat+1] = DATA_deathrate_nonsmoker_educ5[ia][ieduccat]; 
			}
					
		  }

		} // ia+INIT_age<= 50
	}
	infile.close();
	
	
	// Parameter estimates from the HRS analysis: 
	for (int ic=0; ic<13; ic++){ Param_survival_new[ic] = 0.0;  }
	#ifndef MODEL_probs_educ5
			double _param_survival_new[] = { -.0896032,  -1.980246,  .8322108,   2.013633,   .017607,  -.0069537,  -.0176033,  9.738551}; 
			for (int ic=0; ic<8; ic++){ Param_survival_new[ic] = _param_survival_new[ic];  }
	#else 
			double _param_survival_new[] = {  -.0892949,  -1.857563,  .5956055,  .9804853,   1.825659, 2.047653,   .0171861, -.0045391,  -.0082538,  -.0169456,  -.0176481,  -.3042212,  9.939084}; 
			for (int ic=0; ic<13; ic++){ Param_survival_new[ic] = _param_survival_new[ic];  }
	#endif 
	
	
	
	{
	FILE * fdeath ; 
	fdeath  = fopen ("out_data_new/death_hazard.txt", "w"); 
	for(int ia=0; ia < NUM_age90-1; ia++){		
		fprintf(fdeath, "\n age %3d, %f, ", ia+INIT_age, - log( 1.0 - DATA_deathrate_all[ia]) ); 
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			fprintf(fdeath, "%f, ", - log( 1.0 - DATA_deathrate[ia][ieduccat] )); 
		}	
	}
	fclose(fdeath); 
	fdeath  = fopen ("out_data_new/death_rate.txt", "w"); 
	for(int ia=0; ia < NUM_age90-1; ia++){		
		fprintf(fdeath, "\n age %3d, %f, ", ia+INIT_age, DATA_deathrate_all[ia] ); 
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			fprintf(fdeath, "%f, ", DATA_deathrate[ia][ieduccat]); 
		}	
	}
	fclose(fdeath); 	
	fdeath  = fopen ("out_data_new/survival_rate.txt", "w"); 
	for(int ia=0; ia < NUM_age90-1; ia++){		
		fprintf(fdeath, "\n age %3d, %f, ", ia+INIT_age, 1.0-DATA_deathrate_all[ia] ); 
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			fprintf(fdeath, "%f, ", 1.0-DATA_deathrate[ia][ieduccat]); 
		}	
	}
	fclose(fdeath); 	
//	exit(1); 
	}
	
	
	FILE * fdeath ; FILE * felife ; 


	{
	const int educ_cat[] = {8, 12, 13, 16, 17}; 
	int dq=0; double qstock = 0.0; 
	
	dq=0; qstock = 0.0; 
	fdeath  = fopen ("out_data_new/survival_rate_using_dq0_qstock0.txt", "w"); 
	felife  = fopen ("out_data_new/staticLifeExpectancy_using_dq0_qstock0.txt", "w"); 
	for(int ia=0; ia <= NUM_ageMAXT-1; ia++){		
		fprintf(fdeath, "\n age %3d, ", ia+INIT_age ); 
		fprintf(felife, "\n age %3d, ", ia+INIT_age ); 		
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			
			int educ = educ_cat[ieduccat]; int de = 0; 
			double probs = get_survivalProb(ia, educ, qstock, dq, de); double elife = get_lifeExpectancy_static(ia, educ, qstock, dq, de); 
			int current_neduc = 1 + (MAX_initeduc+ia >= 12) + (MAX_initeduc+ia >= 13) + (MAX_initeduc+ia >= 16) + (MAX_initeduc+ia >= 17);
			if ( ieduccat >=  current_neduc ){ probs = -4; elife = -4; }
			fprintf(fdeath, "%f, ", probs); 
			fprintf(felife, "%f, ", elife); 
		}	
	}
	fclose(fdeath); fclose(felife); 

	dq=1; qstock = 1.0; 
	fdeath  = fopen ("out_data_new/survival_rate_using_dq1_qstock1.txt", "w"); 
	felife  = fopen ("out_data_new/staticLifeExpectancy_using_dq1_qstock1.txt", "w"); 
	for(int ia=0; ia <= NUM_ageMAXT-1; ia++){		
		fprintf(fdeath, "\n age %3d, ", ia+INIT_age ); 
		fprintf(felife, "\n age %3d, ", ia+INIT_age ); 		
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			
			int educ = educ_cat[ieduccat]; int de = 0; 
			double probs = get_survivalProb(ia, educ, qstock, dq, de); double elife = get_lifeExpectancy_static(ia, educ, qstock, dq, de); 
			int current_neduc = 1 + (MAX_initeduc+ia >= 12) + (MAX_initeduc+ia >= 13) + (MAX_initeduc+ia >= 16) + (MAX_initeduc+ia >= 17);
			if ( ieduccat >=  current_neduc ){ probs = -4; elife = -4; }
			fprintf(fdeath, "%f, ", probs); 
			fprintf(felife, "%f, ", elife); 
		}	
	}
	fclose(fdeath); fclose(felife); 
	
	dq=0; qstock = 1.0; 
	fdeath  = fopen ("out_data_new/survival_rate_using_dq0_qstock1.txt", "w"); 
	felife  = fopen ("out_data_new/staticLifeExpectancy_using_dq0_qstock1.txt", "w"); 
	for(int ia=0; ia <= NUM_ageMAXT-1; ia++){		
		fprintf(fdeath, "\n age %3d, ", ia+INIT_age ); 
		fprintf(felife, "\n age %3d, ", ia+INIT_age ); 		
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			
			int educ = educ_cat[ieduccat]; int de = 0; 
			double probs = get_survivalProb(ia, educ, qstock, dq, de); double elife = get_lifeExpectancy_static(ia, educ, qstock, dq, de); 
			int current_neduc = 1 + (MAX_initeduc+ia >= 12) + (MAX_initeduc+ia >= 13) + (MAX_initeduc+ia >= 16) + (MAX_initeduc+ia >= 17);
			if ( ieduccat >=  current_neduc ){ probs = -4; elife = -4; }
			fprintf(fdeath, "%f, ", probs); 
			fprintf(felife, "%f, ", elife); 
		}	
	}
	fclose(fdeath); fclose(felife); 
	
	dq=1; qstock = 0.0; 
	fdeath  = fopen ("out_data_new/survival_rate_using_dq1_qstock0.txt", "w"); 
	felife  = fopen ("out_data_new/staticLifeExpectancy_using_dq1_qstock0.txt", "w"); 
	for(int ia=0; ia <= NUM_ageMAXT-1; ia++){		
		fprintf(fdeath, "\n age %3d, ", ia+INIT_age ); 
		fprintf(felife, "\n age %3d, ", ia+INIT_age ); 
		for (int ieduccat =0; ieduccat<neduc_available; ieduccat++){
			
			int educ = educ_cat[ieduccat]; int de = 0; 
			double probs = get_survivalProb(ia, educ, qstock, dq, de); double elife = get_lifeExpectancy_static(ia, educ, qstock, dq, de); 
			int current_neduc = 1 + (MAX_initeduc+ia >= 12) + (MAX_initeduc+ia >= 13) + (MAX_initeduc+ia >= 16) + (MAX_initeduc+ia >= 17);
			if ( ieduccat >=  current_neduc ){ probs = -4; elife = -4; }
			fprintf(fdeath, "%f, ", probs); 
			fprintf(felife, "%f, ", elife); 
		}	
	}
	fclose(fdeath); fclose(felife); 
	
	}
	
	
	
	for(int ia=0; ia < 100; ia++){
	for (int ieduccat =0; ieduccat<NUM_educcat5; ieduccat++){
		DATA_famsize[ia][ieduccat] = 1.0; 
	}}
	

	n = NUM_ismm; 
	
	// initialization, important to initalize off diagnal elements
	for (int i=0; i < n; i++){
	
		TARGET_moments[i] = 0.0;
		TARGET_wgtmoments[i] = 0.0;
	}	
			
	infile.open("out_data_new/smm_moments_NLSY.txt", ifstream::in); 
    if(!infile) {cerr<<"\n Error reading the file smm_moments_NLSY.txt \n"; exit(1); }
    icnt = 0; 
	for(int i=0; i < n && !infile.eof(); i++){
		infile >> TARGET_moments[i];
		icnt ++; 
	}
	infile.close();
	if (icnt!=n) printf("\n Errors in reading smm_moments_NLSY.txt: icnt = %d ", icnt);
	

	infile.open("out_data_new/smm_se_NLSY.txt", ifstream::in); 
    if(!infile) {cerr<<"\n Error reading the file smm_se_NLSY.txt \n"; exit(1); }
    icnt = 0; 
	for(int i=0; i < n && !infile.eof(); i++){
		double temp_se = 0; 
		
		infile >> temp_se; 
		
		TARGET_wgtmoments[i] = (1.0/temp_se ) / temp_se;
		
		icnt ++; 
	}
	infile.close();
	if (icnt!=n) printf("\n Errors in reading smm_se_NLSY.txt: icnt = %d ", icnt);

}




void init_sim(){
// This function generate global variables for value function approximation and equilibrium prices calculation
// All the variables generated here does not depend on parameter values
		
		const double _pi = 3.14159265358979323846264338327950;
		
	#if neps_health==2
		double eps[] = {0.7071067811865475, -0.7071067811865475}; 
    	double weight[] = { 0.8862269254527580, 0.8862269254527580};
    #elif neps_health==3
        double eps[] = {1.224744871391589, 0,  -1.224744871391589};
        double weight[] = { 0.2954089751509193, 1.181635900603677, 0.2954089751509193};
	#elif neps_health==4
		double eps[] = {1.650680123885785, 0.5246476232752903, -0.5246476232752903, -1.650680123885785};
		double weight[] = {0.08131283544724518, 0.8049140900055128, 0.8049140900055128, 0.08131283544724518}; 
    #elif neps_health == 5 
    	double eps[] = {2.020182870456086, 0.9585724646138185, 0.0, -0.9585724646138185, -2.020182870456086};
    	double weight[] = {0.01995324205904591, 0.3936193231522412, 0.9453087204829419, 0.3936193231522412, 0.01995324205904591};
	#elif neps_health == 6 
    	double eps[] = {2.350604973674492, 1.335849074013697, 0.4360774119276165, -0.4360774119276165, -1.335849074013697, -2.350604973674492};
   	 	double weight[] = {0.004530009905508846, 0.1570673203228566, 0.7246295952243925, 0.7246295952243925, 0.1570673203228566, 0.004530009905508846};
   	#elif neps_health == 7
    	double eps[] = {2.651961356835233,1.673551628767471,0.8162878828589647,0.0,-0.8162878828589647,-1.673551628767471,-2.651961356835233};
    	double weight[] = {0.0009717812450995192, 0.05451558281912703, 0.4256072526101278, 0.8102646175568073, 0.4256072526101278, 0.05451558281912703, 0.0009717812450995192}; 
	#elif neps_health == 8
    	double eps[] = {2.930637420257244, 1.981656756695843, 1.157193712446780, 0.3811869902073221, -0.3811869902073221, -1.157193712446780, -1.981656756695843, -2.930637420257244};
    	double weight[] = {0.0001996040722113676, 0.01707798300741348, 0.2078023258148919, 0.6611470125582413, 0.6611470125582413, 0.2078023258148919, 0.01707798300741348, 0.0001996040722113676}; 
	#elif neps_health ==  9
    	double eps[] = {3.190993201781528, 2.266580584531843, 1.468553289216668, 0.7235510187528376, 0.0, -0.7235510187528376, -1.468553289216668, -2.266580584531843, -3.190993201781528};
    	double weight[] = {0.00003960697726326438, 0.004943624275536947, 0.08847452739437657, 0.4326515590025558, 0.7202352156060510, 0.4326515590025558, 0.08847452739437657, 0.004943624275536947, 0.00003960697726326438};    
	#elif neps_health == 10
    	double eps[] = {3.436159118837738, 2.532731674232790, 1.756683649299882, 1.036610829789514, 0.3429013272237046, -0.3429013272237046, -1.036610829789514, -1.756683649299882, -2.532731674232790, -3.436159118837738};
    	double weight[] = {7.640432855232621e-06, 0.001343645746781233, 0.03387439445548106, 0.2401386110823147, 0.6108626337353258, 0.6108626337353258, 0.2401386110823147, 0.03387439445548106, 0.001343645746781233, 7.640432855232621e-06};
	#else 
		printf("\n neps_health invalid, exiting... \n");
		exit(1); 
	#endif 
		
    for (int i=0; i<neps_health; i++){ 		xeps_health[i] = sqrt(2) * eps[i];     weps_health[i] = weight[i]/sqrt(_pi); 	 	}

	/*
    double adj_weps_health[neps_health];
	for (int ih =0; ih < neps_health; ih ++) adj_weps_health [ih] = weps_health[ih];
*/
		
						
		// initialize seeds for shocks
		gsl_rng * rng; 	gsl_rng_env_setup();
		rng = gsl_rng_alloc (gsl_rng_taus);
		//	const gsl_rng_type * T; T = gsl_rng_default;	rng = gsl_rng_alloc (T);	
		
		// shocks in emax function MC integration
		{	
			int seed = 123;		gsl_rng_set( rng, seed ) ; 

			for (int issp = 0; issp < NUM_simsp; issp++){
				for (int i=0; i<NUM_eps; i++){
					simp_rng_eps[issp][i] = gsl_ran_ugaussian(rng); 
					
				}
				
				simp_rng_eps[issp][6] = gsl_rng_uniform(rng); 
				
			}// issp
		}		
		
		// shocks in model simulation. 
		{
			int seed = 123;		gsl_rng_set( rng, seed ) ; 
		
			// Shocks for generating initial distribution of artificial panel			
			for (int ia = 0; ia < NUM_age; ia++){
			for (int iper = 0; iper < NUM_smmPerson; iper++){
				
				for (int i=0; i<NUM_eps; i++){
					sm_rng_eps[ia][iper][i] = gsl_ran_ugaussian(rng); 
				}
				
				sm_rng_eps[ia][iper][6] = gsl_rng_uniform(rng);
				
				
			}
			}
			
		}
		
		
		// Shocks for generating initial distribution of artificial panel 
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
 	 		
 	 		int seed = 12300+100*iunobs;	gsl_rng_set( rng, seed ) ; 
 	 		
 	 		for (int iper = 0; iper < NUM_smmPerson; iper++){
 	 			sm_rng_eps_unobs[iper][iunobs] = gsl_ran_ugaussian(rng); 
				
				if (sm_rng_eps_unobs[iper][iunobs] >  max_ci ) sm_rng_eps_unobs[iper][iunobs] =  max_ci; 
				if (sm_rng_eps_unobs[iper][iunobs] < -max_ci ) sm_rng_eps_unobs[iper][iunobs] = -max_ci; 
			
			}

			// ... only draw the measurement error for health status once for each individual, otherwise the health dynamic equation estimation is contaminated. 		
			for (int iper = 0; iper < NUM_smmPerson; iper++)
                for (int ia = 0; ia < NUM_age+1; ia++){
				sm_rng_eps_measure[ia][iper][iunobs] = gsl_ran_ugaussian(rng);
                }
		}
		
		gsl_rng_free (rng);
}





void init_sm0(int miniper, int maxiper){
// distribution at INIT_age: 


    ifstream infile;   
    
    int n = NUM_DATA_init; 	   
    
    
	double unobs_theta[2][n]; 

//  read data on predicted unobservables from loglike_measures.cpp 	
	infile.open("out_data_new/predicted_thetas.txt"); 	
    if(!infile) {cerr<<"\n Error reading the file loglike_predict_unobs \n"; exit(1); }
    for (int i=0; i < n; i++){	
    	int idnew; 
		infile >> idnew >> unobs_theta[0][i] >> unobs_theta[1][i] ; 
	}
	infile.close(); 
	
		

	double avecostdq[NUM_age+1][n]; double avetaxperpack[NUM_age+1][n]; 
	infile.open("out_data_new/init_distribution_averagecostdq_age15-33.txt"); 	
	if(!infile) {cerr<<"\n\n Restricted access data file (init_distribution_averagecostdq_age15-33.txt) is not available! \n\n"; exit(1); }
    for (int i=0; i < n; i++){	
    	int id; 
    	for (int it=15-INIT_age; it<= 33-INIT_age; it++){
		infile >> id >> avecostdq[it][i] >> avetaxperpack[it][i]; 
		}

		for (int it=33-INIT_age+1; it<NUM_age+1; it++){
		avecostdq[it][i]  = avecostdq[33-INIT_age][i]; 
		avetaxperpack[it][i] = avetaxperpack[33-INIT_age][i]; 
		}
		
		for (int it=15-INIT_age; it<NUM_age+1; it++){
		DATA_avecostdq[it][i]     = avecostdq[it][i]; 
 		DATA_avetaxperpack[it][i] = avetaxperpack[it][i]; 
 		}	
	}
	infile.close(); 
	
	

	
	
	// initialize seeds for random variable draws
	gsl_rng * rng; 	gsl_rng_env_setup();
	rng = gsl_rng_alloc (gsl_rng_taus);
	int	seed = 123;		gsl_rng_set( rng, seed ) ; 
	
		
	int ifirm = 0; 		
	for (int iper = miniper; iper <= maxiper; iper++){	
		
		int ia = 0; 
			
		sm_kstock[ia][iper]       = 0.0; 
		sm_expr[ia][iper]       = 0.0; 		
		
		
		double theta[2]; 

		if (maxiper-miniper+1<NUM_DATA_init){
		
			ifirm = gsl_rng_uniform_int (rng, NUM_DATA_init); 
			
			theta[0] = unobs_theta[0][ifirm]; 
			theta[1] = unobs_theta[1][ifirm];
			
			sm_nlsy_num[iper] = 1; 
			
		} else {
			
			int n_id = floor( (iper-miniper) / (1.0*NUM_DATA_init) ); 

			ifirm = (iper - miniper) - n_id * NUM_DATA_init; 
			
									
			// initialize to posterior mean 
			theta[0] = unobs_theta[0][ifirm]; 
			theta[1] = unobs_theta[1][ifirm];
			
			sm_nlsy_num[iper] = n_id+1; 
			
		}
		
		sm_nlsyid [iper] = DATA_nlsyid[ifirm]; 
		
				
		sm_measure_thetac[iper] = unobs_theta[0][ifirm];  
		
		sm_measure_thetan[iper] = unobs_theta[1][ifirm]; 
		
        sm_isavingspr[iper]  =  1;  
        sm_wgt[iper] = 1.0;
        
		#ifdef DEBUG
		if (ifirm <0 || ifirm >= NUM_DATA_init) {
			printf("\n\n Errors in init_sm0(), it %d \n ... exiting...", ifirm); 	
			exit(1); 
		}
		#endif 
		
		
				
		sm_init_dq[iper][0]        = 0;  
		sm_init_dq[iper][1]        = 0;  
		
		sm_dq_prev[ia][iper]       = 0;   
		
		sm_add   [ia][iper]        = 0; 
		sm_add0 [ia][iper]   	   = (sm_add   [ia][iper] <= 0.01); 
				
		sm_qstock[ia][iper]        = 0;  
		
		
		int hgc = DATA_init_hgc[ifirm];
		if ( hgc > MAX_initeduc ){ 
			hgc = MAX_initeduc; 
		} else if ( hgc < MIN_educ ){
			hgc = MIN_educ; 
		}
				
		sm_educ [ia][iper]       = hgc; 
        sm_de_prev[ia][iper] 	 = DATA_init_de_prev[ifirm];   
		sm_savings[ia][iper]     = DATA_init_net_worth[ifirm]; 	
		

		sm_ieduccatpr[iper]      =  (DATA_init_parents_hgc[ifirm] >= 16);

		
		double min_thetavec[] = {MIN_thetac, MIN_thetan}; double max_thetavec[] = {MAX_thetac, MAX_thetan}; 
		//... bound theta, which depends on parameter value 
		for (int itheta=0; itheta<DIM_theta; itheta++){
		
			if ( theta[itheta] < min_thetavec[itheta] ){ 
				theta[itheta] = min_thetavec[itheta] ; 
			} else if ( theta[itheta] > max_thetavec[itheta] ){ 
				theta[itheta] = max_thetavec[itheta]; 
			}
			
			sm_theta [itheta][iper]   = theta[itheta]; 
		}
				
		for (int ia = 0; ia < NUM_age+1; ia++){
            
            double cprice = DATA_avecostdq[ia][ifirm];
            
            double logcprice = log(cprice);
            if ( logcprice < MIN_logcprice ){
                logcprice = MIN_logcprice;
            } else if ( logcprice > MAX_logcprice ){
                logcprice = MAX_logcprice;
            }
            sm_cprice[ia][iper] = exp(logcprice);
            
			
			// excise tax per pack 
			sm_tax_cprice[ia][iper] = DATA_avetaxperpack[ia][ifirm];  // excise tax per pack
            
        }
		
	
	}

	
	for (int ia = 0; ia < NUM_age-1; ia++){	
	for (int iper = miniper; iper <= maxiper; iper++){
		
		// initialization
		sm_de_prev[ia+1][iper]    = sm_de[ia][iper]; 
		sm_dq_prev[ia+1][iper]    = sm_dq[ia][iper]; 
														
	}}
    
	gsl_rng_free (rng);


}

double get_qnext_itype(double qstock, int dq, int ia, int n){
    
    double qnext = (1.0 - cu_qq) * qstock + 1.0 * dq; 
    if (qnext < 0.0) qnext = 0.0;
    
    
    if (qnext > MAX_qstock[ia+1] && ia < NUM_age ) {
    	qnext = MAX_qstock[ia+1]; 
    }

    return qnext;
}



double get_kstock(int educ, int ia){

	double kstock = ia + INIT_age - 6.0 - educ;
    if (educ < 12) kstock = ia + INIT_age - 6.0 - 12.0;  // start accumulating work expr at age 18.
    if (kstock < 0.0) kstock = 0.0; 
	
	return kstock; 
}


double get_knext(const int ia, const int educ, const double kstock0, const double dk, const int de){
	
	int ia_next = ia+1; 
	int educ_next = educ + de; 
	
	double kstock_next = get_kstock(educ_next, ia_next) ; 
	
    return kstock_next;
}


inline double get_rng_epsw(int educ, double rng_eps){
	
	int ieduc5 = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16); 


	double eps_logwf       = gsl_cdf_gamma_Pinv(rng_eps, sigeps_logwshape[ieduc5], sigeps_logwscale[ieduc5]) + MIN_epsw[ieduc5]; 
	
	return eps_logwf ; 
	
}

inline double get_avg_epsw(int educ){
	
	int ieduc5 = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16); 
	
	return ( sigeps_logwscale[ieduc5]*sigeps_logwshape[ieduc5] + MIN_epsw[ieduc5] );  

}
double get_skill(int educ, double theta[2], int ia, double epsw0){
// ... sigeps_logwscale

    double epsw=epsw0;
  	// ... theoretical mean of the Gamma distribution, with sigeps_logwshape=a, sigeps_logwscale =b 
	double avg_epsw = get_avg_epsw(educ); 
		
	int ieduc4 = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ>=16); 
        
	#ifdef DEBUG 
    if (epsw0 < MIN_epsw[ieduc4]){ printf("\n\n Errors in get_skill: epsw %f < MIN_epsw[ieduc4] %f \n\n ", epsw, MIN_epsw[ieduc4]); exit(1); }
    #endif
	
	double kstock = get_kstock(educ, ia); 
	
	    
	double logskill = clogw_skill[0] * (ia+INIT_age>=18 && ia+INIT_age<=20)*(20.0-(ia+INIT_age))
			+ clogw_skill[1] * kstock 
			+ clogw_skill[2] * kstock*kstock / 100.0 
			+ clogw_skill[3] * (educ - 12) * (educ >= 12)
			+ clogw_skill[4] * (educ - 12) * (educ < 12)
			+ clogw_skill[5] * (educ >= 12 && educ < 16) 
			+ clogw_skill[6] * (educ >= 16) 
			+ clogw_beta[0] 
			+ clogw_beta[1] * kstock * (educ >= 16)   + clogw_ksquare_graduate * kstock*kstock / 100.0 * (educ >= 16) 
			+ clogw_beta[3] * (ia+INIT_age<18)*(18.0-(ia+INIT_age))
			+ theta[0] * clogw_alpha[0][ieduc4] 
			+ theta[1] * clogw_alpha[1][ieduc4] 
			+ clogw_skill_graduate * (educ - 16) * (educ>16); 
			
	
	double skill = exp(logskill) * epsw/avg_epsw; 
    
	    
    #ifdef DEBUG
    if (skill <= 0.0 || isnan(skill) || kstock > MAX_kstock) {
    	printf("\n Errors in get_skill(kstock %.2f, educ %2d, theta, ia %2d, epsw %f) %f < 0.0 || NaN || kstock %f > MAX_kstock %f \n existing... ", kstock, educ, ia, epsw, skill, kstock, MAX_kstock); 
    	exit(1);  
    }
    #endif 
    
	return skill ; 
}


double get_wagerate(const double dk, const int de, const struct_infoset &infoset){

				
	double wagerate = 0.0; 

	if (dk > 0.6) { 
		
		wagerate = infoset.skill;   
	
	} else if (dk > 0.1 && dk < 0.6) {

		wagerate = infoset.skill * exp( clogw_beta[2] * de ) ; 
		
	} 
	

	return wagerate; 
	
}



double get_debtlimit_ptotal(double theta[2], const int ia,  const int educ, const int de, const int dq, const double qstock0){
	
	
	double debt_p = 0.0; 

    if (ia+INIT_age>=18) {

	debt_p = 0.0; 

	int educ_next = educ + de;     
    int ieduc5 = 1*(educ_next==12) + 2 * (educ_next >=13) * (educ_next<=15) + 3 * (educ_next==16) + 4 * (educ_next>16);
	double epsw = MIN_epsw[ieduc5]; 		
	
	// survivalProb depend on qstock 
	double qstock = qstock0; 
	double probs = get_survivalProb(ia, educ, qstock, dq, de); 

				
	for (int ia_next = ia + 1; ia_next < NUM_age; ia_next++){
		
		double skill_next  = get_skill(educ_next, theta, ia_next, epsw);
		
		double inc_wages = skill_next * HOURS; 
		
		double inc_s     = 0.0;  
		
		double incometax = get_incomeTax(inc_wages, inc_s); 
		
		
        #ifdef MODEL_debt0
		
		double cash = ( inc_wages - incometax )  ;     
        
        #else 
        
		double cash = ( inc_wages - incometax - cu_min)  ;  
		
		if (cash <0.0) cash = 0.0; 
		
		#endif 
		
		
		debt_p += ( probs * cash / pow (1.0+R_borrow, ia_next - ia ) ); 
//		debt_p += ( cash / pow (1.0+R_borrow, ia_next - ia ) ); 
			
		
     	double qstock_next  = get_qnext_itype(qstock, dq, ia, 0);
		
		// update survival probability assuming dq_next = dq, de_next = 0 
		probs = get_survivalProb(ia_next, educ_next, qstock_next, dq, 0); 
		
		qstock = qstock_next; 
		
	}
	
	if (debt_p > -1.0* MIN_savings) debt_p = -1.0* MIN_savings; 
	
    }
    
	#ifdef DEBUG 
	if (debt_p < 0.0) printf("\n\n Errors in get_debtlimit_ptotal(%2d, [%.2f, %.2f], %2d), debtlimit = %f < 0.0", educ, theta[0], theta[1], ia, debt_p); 
	#endif 
	
	return debt_p ; 
	
}


double get_debtlimit_endogenous_maxcoef(const int educ_next1, double theta[2], const int ia){

	double debt_p = 0.0;
	
    double thetamax [] = {MAX_thetac, MAX_thetan};
	
	debt_p = get_debtlimit_ptotal(thetamax, ia, educ_next1, 0, 0, 0.0); 
	
	return debt_p ;
		
}





double get_minS(int de, double dk, const double debtlimit_p, const struct_infoset &infoset){
//--- This function returns a lower bound of savings depends on the individual's state
	
	int    educ  = infoset.educ; 
	
	double debt_g = 0.0; 
	
	int itc = educ + de - 13; 
	if (de==1 && itc >=0) {

        
        double debt_before = max( -infoset.savings, 0.0);
        
        debt_g = debt_before + min( max_debtg_e[itc], cbc_tc[itc] - cbc_grants + cbc_roomboard[itc] );

		
		// ... maximum student loan: 35000 (Lochner and Monge-Neranjo 2011)
		double max_debtg1 = max_debtg[0]+(educ+de>16)*max_debtg[1]; 
		if (debt_g > max_debtg1 ) debt_g = max_debtg1; 
		
	}
		
	
	if (de==0 && educ>=16 && infoset.ia + INIT_age <= 25 ) {
		
		debt_g += maxpost_debtg; 
		
	} 
	
	
	#ifdef DEBUG
	if (itc >=8 ) { printf("\n\n Errors in tuition cost and debt ... exiting ... \n\n"); exit(1); }
	#endif 
	
		    
	double min_s = min( -1.0 * debt_g, -1.0 * debtlimit_p); 

		
	#ifdef DEBUG 
		if ( debtlimit_p < -1.0e-4 || min_s > 1.0e-4 ) { printf("\n\n Errors in get_minS: %f , debt_p %f \n \n exiting \n\n", min_s, debtlimit_p ); exit(1); } 
	#endif 
	
	return ( min_s ); 
	
}


double get_ptransfer(int de, double dk, const struct_infoset &infoset){

		
	int age = infoset.ia + INIT_age; 
	
	double ptr = 0.0; 
	
	if ( infoset.ia < NUM_age_ptr )
    {
			
			double logptrstar; 

        logptrstar =  2.245706 * de * (de+infoset.educ >= 13 && de+infoset.educ <= 14)
                + 3.124513 * de * (de+infoset.educ >= 15 && de+infoset.educ <= 16)
                + 1.489327 * de * (de+infoset.educ > 16)
                 -.5464516 * de * (dk>0.1)
                +   .8618642 * (infoset.ieduccatpr==1)
        +  .0456503 * infoset.theta[0]
        +  .1420699 * infoset.theta[1]
        +  .2507142 * infoset.theta[0] * (infoset.ieduccatpr==1) 
        -.2948171 * infoset.theta[1] * (infoset.ieduccatpr==1) 
                 -.2133261 * age
                 -1.92682 * (age==15)
                 -1.759135 * (age==16)
                 -.9232647 * (age==17)
                 +  .019633 * (age==18)
                 -.0562363 * (age==19)
                +  .1332242 * (age>=24)
               -1.810203 * (age>=24) * de * (de+infoset.educ > 13 && de+infoset.educ <= 16)
                +   1.160731 * infoset.de_prev
                +   5.980572
                + infoset.eps_ptr *  2.832729;
							
			ptr = max( exp(logptrstar) - 1.0, 0.0); 

			if (ptr > 100000.0) ptr = 100000.0;  // the maximum from the data
	
	}
	
	return ptr; 
	
}



double get_incomeTax(double inc_wages, double inc_s0){
	
    double tax = TAXrate_y * inc_wages;
	
	return tax; 		

}



inline double getc_subsidy(int ia, int educ, int de){
    

    double c_subsidy = 1.0 * (ia + INIT_age < 18) * c_chi;   
    if (de ==1 && educ < 12) c_subsidy = c_chi;
    
    return (c_subsidy);
}


inline double getc_minexpenditure(int ia, int educ, int de, double c_subsidy){
    
    double roomboard = 0.0;
    if (de == 1 && educ+de>=13) {
        int itc = educ+de - 13; roomboard = cbc_roomboard[itc];
    }
    
    #ifdef consumpEquiv_adj
    
    int ieduccat = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16);
    
    double c_minexpenditure = cu_min*DATA_famsize[ia][ieduccat];
    if ( c_subsidy > 0.0 ) c_minexpenditure = max(0.0, cu_min*DATA_famsize[ia][ieduccat]- c_subsidy);
    
    #else
    
    double c_minexpenditure = cu_min;
    if ( c_subsidy > 0.0 ) c_minexpenditure = max(0.0, cu_min- c_subsidy);
    
    #endif
    
    if ( c_minexpenditure <= roomboard ) c_minexpenditure = roomboard ;
    
    return c_minexpenditure;
}


double get_survivalProb(int ia, int educ0, double qstock, int dq, int de){
	
	
	int educ = educ0 + de;  
	
	
	double probs =  0.0; 
	
	
	const int ageHRS = 51; 
	
	if ( ia + INIT_age < ageHRS ) {
	
	
		int ieduc3 = (educ>=12) + (educ>=16); 
	 	
		probs = 1.0 - DATA_deathrate_nonsmoker_educ3[ia][ieduc3]; 
	

		if (dq >= 1 ) {
		
			probs = 1.0 - DATA_deathrate_smoker_educ3[ia][ieduc3]; 
		
		} 

			
	   #ifdef MODEL_probs_educ5
		
 	   int ieduccat = 0 + 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16);
	
		probs = 1.0 - DATA_deathrate_nonsmoker_educ5[ia][ieduccat]; 
	

		if (dq >= 1 ) {
		
			probs = 1.0 - DATA_deathrate_smoker_educ5[ia][ieduccat]; 
		
		} 
			

		#endif 
		
	} 
	else if (ia + INIT_age >= ageHRS && ia < NUM_ageMAXT-1 ){


		
		#ifndef MODEL_probs_educ5
		
			int ieduccat = 0 + 1*(educ>=12) + 1 * (educ >=16);
			


			double f = Param_survival_new[0] * (ia+INIT_age) 
					+ Param_survival_new[1] * dq 
					+ Param_survival_new[2] * (ieduccat==1)
					+ Param_survival_new[3] * (ieduccat==2) 
					+ Param_survival_new[4] * dq * (ia+INIT_age) 
					+ Param_survival_new[5] * (ieduccat==1) * (ia+INIT_age) 
					+ Param_survival_new[6] * (ieduccat==2) * (ia+INIT_age) 
					+ Param_survival_new[7]; 
					
		
		#else 
			
			
	   		int ieduccat = 0 + 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16);
		

		
			double f = Param_survival_new[0] * (ia+INIT_age) 
					+ Param_survival_new[1] * dq 
					+ Param_survival_new[2] * (ieduccat==1)
					+ Param_survival_new[3] * (ieduccat==2) 
					+ Param_survival_new[4] * (ieduccat==3) 
					+ Param_survival_new[5] * (ieduccat==4) 
					+ Param_survival_new[6] * dq * (ia+INIT_age) 
					+ Param_survival_new[7] * (ieduccat==1) * (ia+INIT_age) 
					+ Param_survival_new[8] * (ieduccat==2) * (ia+INIT_age) 
					+ Param_survival_new[9] * (ieduccat==3) * (ia+INIT_age) 
					+ Param_survival_new[10] * (ieduccat==4) * (ia+INIT_age) 				
					+ Param_survival_new[11] * (dq + qstock>0.0)
					+ Param_survival_new[12]; 
		

		#endif 
		
	
					
		probs = 1.0 - 1.0/(1.0 + exp(f)); 			
		
	} else {
		
		// probs = 0 if NUM_ageMAXT-1
		probs = 0.0; 
	}
	
	
	#ifdef DEBUG
		if ( ia > NUM_ageMAXT-1 || probs < 0.0 || probs > 1.0) { 
		printf("\n Errors in get_survivalProb(ia %d, educ %d, q %f, dq %d) %f > 1.0 \n exiting... \n", ia, educ, qstock, dq, probs); exit(1);
		}
	#endif
	
	return probs; 

}



double get_lifeExpectancy_static(const int ia, const int educ, const double qstock, const int dq, const int de){
		
	double elife = 0.0; 
	
	double xa = 1.0; 
	
	for (int ia_new = ia; ia_new <= NUM_ageMAXT-1; ia_new++) {
		
		// get survival prob at age_new
		double probs = get_survivalProb(ia_new, educ, qstock, dq, de); 
		
		if (probs > 1.0 || probs < 0.0 ) printf("\n Errors in get_lifeExpectancy_static %f \t ia_new %d, educ %d, qstock %f, dq %d \n", probs, ia_new, educ, qstock, dq); 
		
		elife += ( (ia_new + INIT_age ) * (1.0 - probs) * xa ); 
				
		xa *= probs; 	
	}
	
	elife += 1; 
	

	return elife; 
}



double get_lifeExpectancy_dynamic(const int ia0, const int iper){
// Pijoanmas-Rios-Rull2014, pp 2081
// life expectancy (i.e. expected age of death)  
	
	double elife = 0.0; 
	
	double xa = 1.0; 

		double qstock; 
		int dq, de; 
		int    educ; 
		
	for (int ia_new = ia0; ia_new <= NUM_ageMAXT-1; ia_new++) {
		
		if (ia_new <= NUM_age-1) {
			
			qstock = sm_qstock[ia_new][iper]; 
			dq     = sm_dq[ia_new][iper]; 
			educ   = sm_educ[ia_new][iper]; 
			de     = sm_de[ia_new][iper]; 
			
		} else {
		
			qstock = sm_qstock[NUM_age-1][iper]; 
			dq     = sm_dq[NUM_age-1][iper]; 
			educ = sm_educ[NUM_age-1][iper]; 	
			de     = sm_de[NUM_age-1][iper]; 				
		}
		
		// get survival prob at age_new
		double probs = get_survivalProb(ia_new, educ, qstock, dq, de); 
		
		if (probs > 1.0 || probs < 0.0 ) printf("\n Errors in get_lifeExpectancy_dynamic %f \t ia_new %d, educ %d, qstock %f, dq %d \n", probs, ia_new, educ, qstock, dq); 
		
		elife += ( (ia_new + INIT_age ) * (1.0 - probs) * xa ); 
				
		xa *= probs; 	
	}
	
	elife += 1; 
	

	return elife; 
}




double get_delta(double theta[], double qstock){
	
	double delta_theta = delta1[0] * theta[0] + delta1[1] * theta[1]; 
	
	double expterm = - delta0 - delta2[0] * qstock + delta_theta; 
	
	if (expterm > -1.0e-4) expterm = -1.0e-4; 
	
	return exp( expterm ) ;   
	
}




const int NumVfi = ceil(NUM_MGrid/10.0);  
const int NUM_SplusGrid = NUM_MGrid - NumVfi; 
const int NUM_SplusGrid_neg = ceil(0.2*NUM_SplusGrid); 
const int NUM_SplusGrid_pos = NUM_SplusGrid - 1 - NUM_SplusGrid_neg; 
double CURV_SplusGrid_neg = 1.2;   
double CURV_SplusGrid_pos = 1.2; 

double Snext_plusGrid[NUM_age][NUM_ieduc][NUM_SplusGrid];  


#define MAXKINKS 50
#define MAXCOM   50


#define NumChoice_ageYoung 6 
#define NumChoice_ageOld 2 

// FOR each ipt in sparse grid, calculate endogenous M, C, and V. 
double SolutionM_ageYoung_EGM[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_coefEmax][NUM_MGrid], 
       SolutionC_ageYoung_EGM[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_coefEmax][NUM_MGrid], 
       SolutionV_ageYoung_EGM[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_coefEmax][NUM_MGrid]; // 

int    SolutionM_ageYoung_NumKinks[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_coefEmax]; 
double SolutionM_ageYoung_Kinks[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_coefEmax][NUM_MGrid]; 

double Common_ageYoung_MGrid[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_MGrid]; 
double Common_ageYoung_optC[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_MGrid][NUM_coefEmax]; 

// For each m in Common_MGrid, interpolation Coef of sparse grid [thetac, thetan, logcprice, q]
// enable interpolation of consumption function off the grid points 
double CoefV_ageYoung_EGM[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_MGrid][NUM_coefEmax]; //


const int NUM_ageNoEduc = NUM_age-NUM_ageEduc; 

// FOR each ipt in sparse grid, calculate endogenous M, C, and V. 
double SolutionM_ageOld_EGM[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_coefEmax][NUM_MGrid], 
       SolutionC_ageOld_EGM[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_coefEmax][NUM_MGrid], 
       SolutionV_ageOld_EGM[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_coefEmax][NUM_MGrid]; 



double Common_ageOld_MGrid[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_MGrid]; 
double Common_ageOld_optC[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_MGrid][NUM_coefEmax]; 

double CoefV_ageOld_EGM[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_MGrid][NUM_coefEmax]; //


double get_utilityEGM(int de, double dk, const double c, int dq, int ia, int de_prev, int ieduc, int ieduccatpr, double theta[], double qstock, double eps_ue, double eps_uq) {

	int educ = ieduc + MIN_educ; 
	int ieduccat = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16); 
	double uc = pow( c * 0.001, 1.0-cu_g )/(1.0-cu_g) - 1.0/(1.0-cu_g);
	
	int age = ia + INIT_age; 
	double uq = cu_qstock * qstock * qstock + dq * (cu_de_dq *de*dq + cu_q[0] + cu_q[1] * qstock + cu_qe_cat * (educ-MIN_educ) + cu_q_alpha[0] * theta[0] + cu_q_alpha[1] * theta[1] - eps_uq * sigeps_uq  );  
	
	double ue = 0.0; double uke = 0.0; 
	
	
	if (de == 1) {
		
		ue = ( cu_e[0]* (educ + de <= 12) 
				+ cu_e[1]* (educ + de >= 13 && educ + de <= 14) 
				+ cu_e[2]* (educ + de >= 15 && educ + de <= 16)
				+ cu_e[3]* 1.0*( educ + de > 16 )
				- cu_e_rcost * (1.0-de_prev) 
				+ cu_e_alpha[0] * theta[0] + cu_e_alpha[1] * theta[1] 
				+ (ieduccatpr==1)*cu_e_pr
				+ cu_e_a * 1.0 * (age > 23) * ( educ + de <= 16) * (age-22.0)
                + clogw_ndk * ( educ + de >= 10 && educ + de <= 12) * (educ + de - 10)
				- eps_ue * sigeps_ue
				+ cu_ke[0] * (educ + de <= 12) * (18.0-age) 
				+ (cu_e_qstock * qstock)
					);
				
		uke = (dk > 0.1) * ( cu_ke[1] );  
		
	}
	
	#ifdef DEBUG
		if (c < cu_min -1.0 ){ 
			printf("\n\n Errors in get_utilityEGM(%d): c = %f < cu_min %f : %f,  \n ... exiting ...", ia, c, cu_min, cu_min*DATA_famsize[ia][ieduccat]); // exit(1);  
		}
		if (uc < 0.0-1.0e-8) {
			printf("\n\n Errors in get_utilityEGM(%d): cu = %.8f \t, c = %f, cu_min = %f : %f \n ... exiting ...", ia, uc, c, cu_min, cu_min*DATA_famsize[ia][ieduccat]); exit(1);  		
		}
		if (de+dk>1.5) {
			printf("\n\n Errors in get_utilityEGM(%d): de %d, dk %f \n ... exiting ...", ia, de, dk); exit(1);  
		}
		if (educ < MIN_educ || educ > MAX_educ ){
			printf("\n\n Errors in get_utilityEGM(%d): educ %d \n ... exiting ...", ia, educ); exit(1);  		}
	#endif	

	
	return ( uc + uq + (ue + uke) ); 
		

}
inline double get_dUdC_EGM(const double c, const int ia, const int educ){

	#ifdef DEBUG 
		
		if (c<1000.0) { printf("\n Errors in get_dUdC_EGM(), c = %f < 1000.0 \n exiting...\n", c); exit(1); }
	
	#endif
		
	
	return ( pow( c * 0.001, -cu_g) ); // * (0.001/DATA_famsize[ia][ieduccat]) );

} 
inline double inv_dUdC_EGM(const double du, const int ia, const int educ){
	
	#ifdef DEBUG 
		
		if (du<1.0e-10) { printf("\n Errors in inv_dUdC_EGM(), du = %f < 1.0e-10 \n exiting...\n", du); exit(1); }
	
	#endif
	

	double c = pow( du, -1.0/cu_g ) / 0.001; // / (0.001/DATA_famsize[ia][ieduccat]) ;
	
	return c; 
} 



double get_PDV_EGM_new(int ia, int ieduc, double theta[], double kstock, double health, double savings0, int itype, const int dq0, double qstock0, double &EdV_65){
		
	#ifdef DEBUG
	
	int T_next = INIT_age + NUM_age - 1 + 1; 
	
	if (ia != NUM_age) { printf("\n Errors in get_PDV: Terminal decision age = %2d  \n exiting \n ", INIT_age + NUM_age - 1); exit(1); }
    if (savings0 <  -1.0e-6){ printf("\n Errors in get_PDV: Terminal savings %f < 0.0 \n exiting \n ", savings0); exit(1); }
 	if (T_next !=65) { printf("\n Errors in get_PDV: Terminal age %d != 65 \n exiting \n ", T_next); exit(1); } 
 	
	#endif 
	
	
	int dq = dq0; 
	

	double vecY_retire[] = { 9125.378 - health * dq,   11200.89 - health * dq,  11325.07 - health * dq, 11232.91 - health * dq, 11590.56 - health * dq  };
	
	
   	int educ = ieduc+MIN_educ; 
	int ieduccat = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16); 	   	
	
	 double r = R_lend*(savings0>0.0) + R_borrow*(savings0<0.0); 

	 const double price_c = 1.0/(1.0+R_lend); 	

	 double num_retire = NUM_ageMAXT-1 - NUM_age; 

	 double pdv_Y      =  vecY_retire[ieduccat] * (1.0 - pow( price_c, num_retire +1)) / (1.0 - price_c); 
	 
		 
	
	 double qstock = qstock0; 
	 double total_coef = 0.0; 
	 double compound_sdelta = 1.0; 
	 double PDV_retire = 0.0; 
	 double savings  = savings0; double c65_retire = 0.0; 
	 
	 double coefvec[NUM_ageMAXT - NUM_age+1]; 
	 
	 total_coef = 0.0; 
	 compound_sdelta = 1.0; qstock = qstock0; dq = dq0; 
	 for (int ia = NUM_age; ia < NUM_ageMAXT; ia++){

		double compound_price_c = pow( price_c, ia-NUM_age); 
		coefvec[ia-NUM_age] =  pow( compound_price_c / compound_sdelta, -1.0/cu_g); 
		total_coef += (compound_price_c * coefvec[ia-NUM_age]) ; 
		
		// update compound_probsdelta: 
		double probs = get_survivalProb(ia, educ, qstock, dq, 0); 
		double delta = get_delta(theta, qstock); 	
	 	compound_sdelta *= (probs * delta); 
			
    	// update qstock 
     	qstock  = get_qnext_itype(qstock, dq0, ia, itype);
		
	 } 	 
	  
	 // marginal value of wealth at optimal (allowing borrowing across periods)
 	 double pdv_wealth = (1.0+r)*savings0 + pdv_Y; 
 	 
		 
	if (pdv_wealth > 0.0 ) {
		
	 //		 double total_coef_check = 0.0; 
	 // ... 
	 PDV_retire = 0.0; 
	 savings  = savings0;  c65_retire = 0.0;  
	 compound_sdelta = 1.0; qstock = qstock0; dq = dq0; 
	 for (int ia = NUM_age; ia < NUM_ageMAXT; ia++){
				
		double  c_retire = 	pdv_wealth  * coefvec[ia-NUM_age] / total_coef; 

	 	double cash = (1.0+r)*savings + vecY_retire[ieduccat]; 
 		if (c_retire < cu_min ) c_retire = cu_min; 
 		if (c_retire > cash ) c_retire = cash; 
		 		
 		double uc_retire = get_utilityEGM(0, 0, c_retire, dq, ia, 0, ieduc, 0, theta, qstock, 0.0, 0.0); 
     	
		// update compound_probsdelta: 
		double probs = get_survivalProb(ia, educ, qstock, dq, 0); 
		double delta = get_delta(theta, qstock); 	
	 	compound_sdelta *= (probs * delta); 
		 
     	PDV_retire += (compound_sdelta*uc_retire);  
		
     	// update qstock and snext 
     	qstock  = get_qnext_itype(qstock, dq0, ia, itype);
     	savings = cash - c_retire; 
     	     	     	
		if (ia==NUM_age) c65_retire = c_retire; 	     	
	 } 	 


	double Edu_65 = get_dUdC_EGM(c65_retire, NUM_age, educ); 
	
	EdV_65 =  (1.0+r)*Edu_65; 
	
	
	} else {
	
		// pdv_wealth < 0 could occur if smoking expenditure is too high compared to savings0 and ssi income 
		
		PDV_retire = - 1000000000.0; 
		
		EdV_65 = 0.0; 
	
	} 
	
	return PDV_retire; 
	
}



void getpar_XWEPS_logwf(){		
		// global variable used for numerical integration over earnings shocks
    
		    VecDoub xeps(NUM_EPS_earnings); 
		    VecDoub weps_earnings(NUM_EPS_earnings); 
			for (int ieduc5 = 0; ieduc5 < 5; ieduc5 ++){
						    	
				double shape = sigeps_logwshape[ieduc5]; 
				double scale = sigeps_logwscale[ieduc5]; 
				
		    	Doub alf_earnings = shape -1.0; 
		    	
		    	// ... Gauss-Laguerre:  0 < x < infty => earnings shock, can be directly used to approximate gamma distribution with scale = 1 
		    	gaulag(xeps, weps_earnings, alf_earnings); 
				
				for (int ie =0; ie < NUM_EPS_earnings; ie++) { 	
		    		XEPS_logwf[ieduc5][ie]     = xeps[ie] * scale + MIN_epsw[ieduc5]; 
		    		WEPS_logwf[ieduc5][ie]     = weps_earnings[ie] * pow( xeps[ie], -alf_earnings) * exp( xeps[ie] ) * gsl_ran_gamma_pdf (xeps[ie], shape, 1.0); 
		    	}
																								
			}

				
}	



void getpar(const gsl_vector* x) {

	int ivec = 0;

	cu_g =  (gsl_vector_get(x, ivec) ) ;  				ivec++;

    clogw_skill_graduate = ( gsl_vector_get(x, ivec) ); 		ivec++;    
	
	delta0 =  (gsl_vector_get(x, ivec) ) ;  				ivec++;
	
    delta1[0] =  (gsl_vector_get(x, ivec) ) ;  				ivec++;
	
    delta1[1] = (gsl_vector_get(x, ivec) ) ;  				ivec++;

    delta2[0] = (gsl_vector_get(x, ivec) ) ;  				ivec++;
	
	
 	cu_e[0] = ( gsl_vector_get(x, ivec ) );								ivec++; 

	cu_e[1] = ( gsl_vector_get(x, ivec ) );			ivec++;
    
    cu_e[2] = ( gsl_vector_get(x, ivec ) );						ivec++;

    cu_e[3] = ( gsl_vector_get(x, ivec ) );						ivec++;


	cu_e_rcost = ( gsl_vector_get(x, ivec) ); 						ivec++;  	
	
    cu_e_alpha[0]    = ( gsl_vector_get(x, ivec ) ); 		ivec++;  

	cu_e_alpha[1]    = ( gsl_vector_get(x, ivec ) ); 		ivec++;  
	
    cu_e_pr   = ( gsl_vector_get(x, ivec) );							ivec++;  	
	
 	cu_e_a   = ( gsl_vector_get(x, ivec) );				ivec++;
	
	
	       
 	sigeps_ue = ( gsl_vector_get(x, ivec ) ); 			ivec++;
	
		
    cu_ke[1] = ( gsl_vector_get(x, ivec) ); 			ivec++;
	
    cu_ke[0] = ( gsl_vector_get(x, ivec) ); 			ivec++;
		

	cu_q[0] = gsl_vector_get(x, ivec ); 				ivec++;  	
	
	cu_q[1] = ( gsl_vector_get(x, ivec) ); 			ivec++;  
	
	cu_qe_cat 		= ( gsl_vector_get(x, ivec) ); 	 	ivec++;  	

	cu_qq = gsl_vector_get(x, ivec) ;  ivec++;  	
	
	cu_q_alpha[0]   = ( gsl_vector_get(x, ivec ) ); 		ivec++;   
	cu_q_alpha[1]   = ( gsl_vector_get(x, ivec ) ); 		ivec++;  

	cu_qstock 		    = ( gsl_vector_get(x, ivec ) ); 	 	ivec++;  	

	sigeps_uq  		= ( gsl_vector_get(x, ivec ) ); 		ivec++;  

	

	
	// human capital and wages
	clogw_skill[0]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
		
	clogw_skill[1]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
	
    clogw_skill[2]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
		
    clogw_skill[3]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
		
 	
    clogw_skill[4]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
	clogw_skill[5]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
	clogw_skill[6]  = ( gsl_vector_get(x, ivec) ); 		ivec++;  	 
		


  	
	// wages

	clogw_beta[0]    =   gsl_vector_get(x, ivec ); 			ivec++;  	

	clogw_beta[1]  =  (  gsl_vector_get(x, ivec) ) ; 		ivec++;  	
	    
    clogw_beta[2]  =  (  gsl_vector_get(x, ivec ) ); 		ivec++;  	
	    
    clogw_beta[3]  =  (  gsl_vector_get(x, ivec ) ); 		ivec++;  	

	clogw_ndk      =  (  gsl_vector_get(x, ivec ) ); 		ivec++;  	
	
	
	clogw_alpha[0][0]   =  ( gsl_vector_get(x, ivec ) ); 		ivec++; 
	clogw_alpha[1][0]   =  ( gsl_vector_get(x, ivec ) ); 		ivec++; 
	
	for (int ieduccat = 1; ieduccat < 4; ieduccat ++) {
	clogw_alpha[0][ieduccat] =  clogw_alpha[0][0];     
	clogw_alpha[1][ieduccat] =  clogw_alpha[1][0];     
 	}
 	
	

	// ... wage shocks 
    sigeps_logwshape[0] = - gsl_vector_get(x, ivec )*10.0;     ivec++;
    sigeps_logwshape[1] = sigeps_logwshape[0];
    sigeps_logwshape[2] = sigeps_logwshape[1];
    
    sigeps_logwshape[3] = sigeps_logwshape[0]; 
    sigeps_logwshape[4] = sigeps_logwshape[3] ;

    
    sigeps_logwscale[0]  = -gsl_vector_get(x, ivec )*100.0; 			ivec++;
    sigeps_logwscale[1]  = sigeps_logwscale[0];
    sigeps_logwscale[2]  = sigeps_logwscale[1];
    
    sigeps_logwscale[3]  = sigeps_logwscale[0]; 
    if (sigeps_logwscale[3] < sigeps_logwscale[0]) sigeps_logwscale[3] = sigeps_logwscale[0]; 
    sigeps_logwscale[4]  = sigeps_logwscale[3]; 
    
		
	cu_e_qstock = ( gsl_vector_get(x, ivec ) ); 		ivec++;  


}



void getvec(gsl_vector* x) {

	int ivec = 0;

    gsl_vector_set(x, ivec, ( cu_g  ) );	ivec++;

    gsl_vector_set(x, ivec, ( clogw_skill_graduate ) ); ivec++;    
    
    gsl_vector_set(x, ivec, ( delta0  ) );	ivec++;
	
    gsl_vector_set(x, ivec, ( delta1[0] ) );	ivec++;

    gsl_vector_set(x, ivec, ( delta1[1] ) );	ivec++;
	
	gsl_vector_set(x, ivec, ( delta2[0] ) );	ivec++;
	

	
 	gsl_vector_set(x, ivec, ( cu_e[0] ) ); 					ivec++;  

	gsl_vector_set(x, ivec, ( cu_e[1] ) ); 				ivec++;  

    gsl_vector_set(x, ivec, ( cu_e[2]  ) ); 					ivec++;

    gsl_vector_set(x, ivec, ( cu_e[3]  ) ); 					ivec++;


	gsl_vector_set(x, ivec, ( cu_e_rcost ) ); 		ivec++;  	
	  
    gsl_vector_set(x, ivec, ( cu_e_alpha[0] ) ); 		    ivec++;  
	
 	gsl_vector_set(x, ivec, ( cu_e_alpha[1] ) ); 		    ivec++;  
	
    gsl_vector_set(x, ivec, (cu_e_pr ) );   			ivec++;  	

 	gsl_vector_set(x, ivec, (cu_e_a ) );   				ivec++;

 		
    gsl_vector_set(x, ivec, (sigeps_ue) ); 		ivec++;


 	gsl_vector_set(x, ivec, ( cu_ke[1] ) ); 			ivec++;  	

 	gsl_vector_set(x, ivec, ( cu_ke[0] ) ); 			ivec++;  	
	
		
	gsl_vector_set(x, ivec, cu_q[0] ); 					ivec++;  	

	gsl_vector_set(x, ivec, (cu_q[1]) ); 		ivec++;  	
	
	gsl_vector_set(x, ivec, ( cu_qe_cat    ) ); 		ivec++;  	

	gsl_vector_set(x, ivec, cu_qq );	ivec++; 	

 	gsl_vector_set(x, ivec, ( cu_q_alpha[0] ) ); 			ivec++;   
 	gsl_vector_set(x, ivec, ( cu_q_alpha[1] ) ); 			ivec++;   

	gsl_vector_set(x, ivec, ( cu_qstock   ) ); 			ivec++;  	
	
 	gsl_vector_set(x, ivec, (sigeps_uq ) ); 		ivec++;  
	

	
	
	// human capital and wages
	gsl_vector_set(x, ivec, ( clogw_skill[0]) ); 		ivec++;  	
	
	gsl_vector_set(x, ivec, ( clogw_skill[1]) ); 		ivec++;  	

    gsl_vector_set(x, ivec, ( clogw_skill[2]) ); 		ivec++;  	
	
	gsl_vector_set(x, ivec, ( clogw_skill[3]) ); 		ivec++;  	

    gsl_vector_set(x, ivec, ( clogw_skill[4]) ); 		ivec++;  	
	gsl_vector_set(x, ivec, ( clogw_skill[5]) ); 		ivec++;  	
	gsl_vector_set(x, ivec, ( clogw_skill[6]) ); 		ivec++;  	

		

	// wages

    gsl_vector_set(x, ivec,      clogw_beta[0] ); 		    ivec++;  	

    gsl_vector_set(x, ivec,  ( clogw_beta[1] ) ); 		    ivec++;  	

	gsl_vector_set(x, ivec,  ( clogw_beta[2] ) ); 		    ivec++;  	

	gsl_vector_set(x, ivec,  ( clogw_beta[3] ) ); 		    ivec++;  	

	gsl_vector_set(x, ivec,  ( clogw_ndk ) ); 		    ivec++;  		
	
	
	gsl_vector_set(x, ivec,    ( clogw_alpha[0][0] ) ); 		ivec++; 
	gsl_vector_set(x, ivec,    ( clogw_alpha[1][0] ) ); 		ivec++; 
	
	
	// wage shocks
    gsl_vector_set(x, ivec, (-sigeps_logwshape[0]/10.0) );     ivec++;
    
    gsl_vector_set(x, ivec, (-sigeps_logwscale[0]/100.0) ); 	ivec++;
	
	gsl_vector_set(x, ivec,    ( cu_e_qstock ) ); 		ivec++;  


}



void get_coefEmaxstate(){
// must obtain after the parameter values
	
	for (int ieduc = 0; ieduc < NUM_ieduc; ieduc++){			

    int educ = ieduc+MIN_educ;
        
        if (educ <= 12){
            MAX_debt_gsl[ieduc] = 0.0;
        } else {
            
            int itc = educ-13;
            double debt_g = 0.0;
            for (int i=0; i<= itc; i++) debt_g += max_debtg_e[i];
            
            MAX_debt_gsl[ieduc] = min( max_debtg[0]+(educ>16)*max_debtg[1], debt_g );
            

        }
	
	if (MAX_debt_gsl[NUM_ieduc-1] > -1.0 * MIN_savings){
		printf("\n\n Error: MAX_debt_gsl[NUM_ieduc-1] %f > - MIN_savings %f \n Must decrease MIN_savings \n\n exiting ... \n ", MAX_debt_gsl[NUM_ieduc-1], MIN_savings); exit(1); 
	}
	
	int ieduc4 = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ>=16); 
    int ieduc5 = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16);
        
    double theta [] = {MAX_thetac, MAX_thetan};
	
	if (clogw_alpha[1][ieduc4] < 0.0) theta [1] = MIN_thetan;  
				
		{
			
		// ... initialization ... 	
		MAX_debt [0][ieduc] = 1.0; 		MAX_asset[0][ieduc] = 1.0; 
		

		MAX_debt[1][ieduc]  = 5000.0; 	MAX_asset[1][ieduc] = 100000.0;
		

		MAX_debt[2][ieduc]  = max(3000.0, MAX_debt_gsl[ieduc]); 	MAX_asset[2][ieduc] = 100000.0;
		
				

		for (int ia = 3; ia < NUM_age; ia++){
			

			double maxdb_private = get_debtlimit_endogenous_maxcoef(educ, theta, ia-1)+100.0;
			
			double maxdb_gsl = MAX_debt_gsl[ieduc]; 
			if (ia >= NUM_ageEduc) maxdb_gsl = ( 1.0-(ia-NUM_ageEduc)/30.0 ) * MAX_debt_gsl[ieduc];  
			
            MAX_debt[ia][ieduc] = max(maxdb_private, maxdb_gsl);
			if (MAX_debt[ia][ieduc] < 3000.0) MAX_debt[ia][ieduc] = 3000.0; 
        			            
			if (ia + INIT_age <= 20 || ieduc > ia + (MAX_initeduc - MIN_educ) ){
				
				MAX_asset[ia][ieduc] = 100000.0; 
				
			} else {
				
				MAX_asset[ia][ieduc]= (1.0+R_lend)*MAX_asset[ia-1][ieduc] + 20000.0; // * kstock;
			}
			
			if (MAX_debt[ia][ieduc] > -MIN_savings) MAX_debt[ia][ieduc] = -MIN_savings;
			if (MAX_asset[ia][ieduc] > MAX_savings[ieduc5]) MAX_asset[ia][ieduc] = MAX_savings[ieduc5];
			
		} 
		


		MAX_debt [NUM_age][ieduc] = MAX_debt [NUM_age-1][ieduc]; 		MAX_asset[NUM_age][ieduc] =MAX_savings[ieduc5];
		
		

		for (int ia = 0; ia < NUM_age; ia++){

			int isg=0; 
			for (int i=NUM_networth_neg; i>=1; i--){
				networthGrid[ia][ieduc][isg] = -MAX_debt[ia][ieduc] * pow( 1.0*i/NUM_networth_neg, curv_networth_neg) + 0.0; isg++; 
			}
			networthGrid[ia][ieduc][isg] = 0.0; isg++;	
			for (int i=1; i<=NUM_networth_pos; i++){
				networthGrid[ia][ieduc][isg] = MAX_asset[ia][ieduc] * pow( 1.0*i/NUM_networth_pos, curv_networth_pos) + 0.0; isg++;
			} 
						
			// specify Splus Grid here: 
			isg=0; 
			for (int i=NUM_SplusGrid_neg; i>=1; i--){
				Snext_plusGrid[ia][ieduc][isg] = -MAX_debt[ia+1][ieduc] * pow( 1.0*i/NUM_SplusGrid_neg, CURV_SplusGrid_neg) + MAX_debt[ia+1][ieduc] + 1.0e-4; isg++; 
			}
			Snext_plusGrid[ia][ieduc][isg] = MAX_debt[ia+1][ieduc]; isg++;	
			for (int i=1; i<=NUM_SplusGrid_pos; i++){
				Snext_plusGrid[ia][ieduc][isg] = MAX_asset[ia+1][ieduc] * pow( 1.0*i/NUM_SplusGrid_pos, CURV_SplusGrid_pos) + MAX_debt[ia+1][ieduc]; isg++;
			} 

			if ( fabs( Snext_plusGrid[ia][ieduc][0] ) > 1.0e-4){ printf("\n Errors in Splus_Grid, please recheck!!! \n exiting \n "); exit(1); } 
				
		}

	} // max_debt, max_asset	
	
	} //ieduc
	
	
	// ... must call first, before calculating MAX_debt an MAX_asset
	MAX_qstock[0] = 0.0001; // the number is positive in order to expand the state space. 
	MAX_qstock[1] = 1.0; 
	
	for (int ia = 2; ia < 100; ia++){
		
		MAX_qstock[ia] = 1.0 + (1.0-cu_qq) * MAX_qstock[ia-1];  // 1.0 * ia; 
	}
	
	
	for (int ia = 0; ia < NUM_age; ia++){	
	for (int ieduc = 0; ieduc < NUM_ieduc; ieduc++){	
        
		MAX_coefEmaxstate[ia][ieduc][0] = exp(MAX_thetac); 
		MIN_coefEmaxstate[ia][ieduc][0] = exp(MIN_thetac); 
		
		MAX_coefEmaxstate[ia][ieduc][1] = exp(MAX_thetan); 
		MIN_coefEmaxstate[ia][ieduc][1] = exp(MIN_thetan); 
		
		
		MIN_coefEmaxstate[ia][ieduc][2] = MIN_logcprice;    		
		MAX_coefEmaxstate[ia][ieduc][2] = MAX_logcprice;
		
		
		// qstock 
		MIN_coefEmaxstate[ia][ieduc][3] = 0.0; 				
		MAX_coefEmaxstate[ia][ieduc][3] = MAX_qstock[ia]; 	
		
        
		#ifdef DEBUG 
		if ( ia >= 0 && ( (MAX_coefEmaxstate[ia][ieduc][0]-MIN_coefEmaxstate[ia][ieduc][0]) <= 0.01 || (MAX_coefEmaxstate[ia][ieduc][1]-MIN_coefEmaxstate[ia][ieduc][1]) <= 0.01  )) printf("\n Errors in init_coefEmax: ia %2d, ieduc %2d, logcprice [%f, %f], thetac [%.2f, %.2f], thetan [%.2f, %.2f] \n ... exiting \n", ia, ieduc, MIN_coefEmaxstate[ia][ieduc][2], MAX_coefEmaxstate[ia][ieduc][2], MIN_coefEmaxstate[ia][ieduc][0],MAX_coefEmaxstate[ia][ieduc][0], MIN_coefEmaxstate[ia][ieduc][1], MAX_coefEmaxstate[ia][ieduc][1]);
		#endif 
	}}	

}


void init_coefEmax(){
				
	//!!! must first call init_par() or getpar() !!!
	get_coefEmaxstate(); 
	
	
	//!!!! Initialization of value function coefficients
	for (int ia = NUM_age-1; ia >= 0; ia--){
	
		for (int de_prev = 0; de_prev < NUM_de; de_prev++){  
		for (int ieduc = 0; ieduc < NUM_ieduc; ieduc++){
			
			for (int isavingspr = 0; isavingspr < NUM_isavingspr; isavingspr++){
			for (int ieduccatpr = 0; ieduccatpr < NUM_ieduccatpr; ieduccatpr++){
			for (int isavings = 0; isavings < NUM_networth; isavings++){
			for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
				
				#ifndef SOLVE_EGM
				coefEmax[ia][de_prev][ieduc][isavingspr][ieduccatpr][isavings][ipt] = 0.0;	
				#endif 
				
				#ifdef MODEL_dEmax
				coefdEmax[ia][de_prev][ieduc][isavingspr][ieduccatpr][isavings][ipt] = 0.0;	
				#endif 
				
				#ifdef MODEL_emaxQ0
				coefEmax_Q0[ia][de_prev][ieduc][isavingspr][ieduccatpr][isavings][ipt] = 0.0;	
				#endif 
			}}
			}}

			
		}}
	}
	#ifndef SOLVE_EGM
	const int n =(NUM_age)*(NUM_de)*(NUM_ieduc)*NUM_isavingspr;
    printf("\n\n the dimension of Emax function for given theta: %d \n\n", n);  
	#endif 
	
	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	// RULE=1 "Clenshaw Curtis Closed Fully Nested Rule"	
	// Number of collocation points on each dimension of state variables
	// max{i_d, i_1+i_2+...+i_d <= ndim+level} = 1+level
	// maximum order of polynomial: 2^{i_d-1}+1 = 2^{level}+1 = 5 in current case
		
	int ndim = DIM_smolyak; 
	
	int level = my_max(DIM_smolyak, LEVEL_smolyak_vector); 
	
	int order = levels_index_size(1, level, RULE_smolyak);		
	
	int npt   = levels_index_size(ndim, level, RULE_smolyak);	// number of distinct collocation points
	
	LEVEL_smolyak = level; 
	ORDER_smolyak = order; 
	
	int min_level = my_min(DIM_smolyak, LEVEL_smolyak_vector); 
	
	#ifdef DEBUG
	printf("\n\n Initializing state space using Sparse Grid Method, ndim = %d, Level = %d \n\n", ndim, level);
	printf(" Clenshaw Curtis Closed Fully Nested Rule \n"); 
	printf(" The product grid has order %d, level %d \n", order, level); 
	printf(" The number of distinct points in a sparse grid %d \n\n", npt);
	#endif	
	

	// FIND UNIQUE SET OF POLYNOMIALS OF TOTAL DEGREE level IN ALL dimENSION 

	#ifdef DEBUG_smolyak
	printf("\n Sparse grid points = %3d, level = %d, max(ndim, level+1) = %d \n", NUM_coefEmax, LEVEL_smolyak, max(ndim, LEVEL_smolyak+1));	
	#endif

	// T_(l_1) * T_(l_2) * ... * T_(l_dim) for P^I

	int l_min = max ( LEVEL_smolyak + 1 - ndim, 0 );
	
	if (min_level == level) {
	
	ChebyIndex = my_allocate_int(npt, ndim);
	ChebyIndex = getChebyIndex(l_min, LEVEL_smolyak, npt, ndim); 
	
	} else {
		
	int **cheby_iso; cheby_iso = my_allocate_int(npt, ndim); 
	cheby_iso = getChebyIndex(l_min, LEVEL_smolyak, npt, ndim);  
	
		
	#ifdef DEBUG_smolyak
	printf("\n ilevel = %d, [I_min = %d , I_max = %d], npt=%d \n", level, l_min+ndim, LEVEL_smolyak+ndim, npt);
	for (int inm =0; inm< npt; inm++){
		printf("\n ipt %3d, ipoly = %3d,\t Cheby  = ", inm, my_sum(ndim, cheby_iso[inm]));
		for (int idim=0; idim< ndim; idim++){
			printf(" %d \t ", cheby_iso[inm][idim] );	
		}
	}
	printf("\n");
	#endif	
	
			
	// ... prepare for Smolyak Anisotrop 
	int vector_good[npt]; int npt_aniso; 
	getChebyIndex_Anisotrop_size(npt, ndim, cheby_iso, LEVEL_smolyak_vector, vector_good, npt_aniso); 
	
	printf("\n npt_aniso %2d \n", npt_aniso);
	 
	int npt0 = npt; 
		
	npt = npt_aniso; 
	
	
	ChebyIndex = my_allocate_int(npt, ndim);
	
//	ChebyIndex = getChebyIndex(l_min, LEVEL_smolyak, npt, ndim);  
		
	ChebyIndex = getChebyIndex_Anisotrop(npt0, ndim, cheby_iso, npt_aniso, vector_good); 
	
	my_free(cheby_iso, npt); 
		
	}
	
	
	#ifdef DEBUG_smolyak
	printf("\n ilevel = %d, [I_min = %d , I_max = %d], npt=%d \n", level, l_min+ndim, LEVEL_smolyak+ndim, npt);
	for (int inm =0; inm< npt; inm++){
		printf("\n ipt %3d, ipoly = %3d,\t Cheby  = ", inm, my_sum(ndim, ChebyIndex[inm]));
		for (int idim=0; idim< ndim; idim++){
			printf(" %d \t ", ChebyIndex[inm][idim] );	
		}
	}
	printf("\n");
	#endif	
		
		
	GRID_smolyak   = my_allocate_db(npt, ndim);
	
	double grid_point[ndim*npt];

	Smolyak_Grid(ndim, level, npt, ChebyIndex, grid_point); 
	
	for (int ip=0; ip<npt; ip++){
		for (int idim =0; idim<ndim; idim++){
		
			GRID_smolyak[ip][idim] = grid_point[idim + ip * ndim];
		}
	}
	
	
	#ifdef DEBUG_smolyak
	printf("\n ilevel = %d, [I_min = %d , I_max = %d], npt=%d \n", level, l_min+ndim, LEVEL_smolyak+ndim, npt);
	for (int inm =0; inm< npt; inm++){
		printf("\n ipt %3d, ipoly = %3d,\t Cheby  = ", inm, my_sum(ndim, ChebyIndex[inm]));
		for (int idim=0; idim< ndim; idim++){
			printf(" %d (%.4f) \t ", ChebyIndex[inm][idim], GRID_smolyak[inm][idim] );	
		}
	}
	printf("\n");
	#endif		
		
	
	#ifdef DEBUG_smolyak
	std:: cout<< "\n Sparse grid is initialized \n\n";
	#endif
		
	
	// ... to avoid repeated inversion of matrix: 
	
	POLYINV_smolyak = gsl_matrix_alloc(npt, npt); 
	
	struct_smolyak smolyak; 
	smolyak.PolyInvert(ChebyIndex, GRID_smolyak, POLYINV_smolyak); 
	

	
	if ( NUM_coefEmax != npt ){
	
		printf("\n\n Error in NUM_coefEmax = %d != npt = %d \n ... exiting \n", NUM_coefEmax, npt);
		exit(1);  
	}
	
	
	#ifdef DEBUG_smolyak
	
	double **emaxstates; 
	emaxstates = my_allocate_db(NUM_coefEmax, DIM_smolyak);
	
	int ia = 1; int ieduc = 12 - MIN_educ; 
	
	// states =  (states_norm + 1.0)(state_max-state_min)/2.0 + state_min
	for (int ipt=0; ipt < NUM_coefEmax; ipt++){
		
		for (int is=0; is<ndim; is++){
			emaxstates[ipt][is] =  (GRID_smolyak[ipt][is] + 1.0)*(MAX_coefEmaxstate[ia][ieduc][is]-MIN_coefEmaxstate[ia][ieduc][is])/2.0 + MIN_coefEmaxstate[ia][ieduc][is];
		}
		
		#ifdef DEBUG_smolyak
			for (int idim=0; idim< ndim; idim++) std::cout<< emaxstates[ipt][idim] << "\t";
			std::cout <<"\n"; 
		#endif
		
	}
			
	printf("\n\n Finishing debug smolyak method \n\n exiting..\n\n"); 
	
	my_free(emaxstates, NUM_coefEmax); 
	
	exit(1); 
	
	#endif 
	
	
}


void init(){

    init_data(); 
	
	init_par(); 
	initpar_loglike(); 
	
	
	init_data(); 

  	ifstream infile; 			
    infile.open("estpar_iter.txt"); 	    readpar(500, infile);	infile.close();

    init_sim(); 
    	    
    init_coefEmax(); 
	
	cout << "\n\n init() completed \n\n"; 	    
		
}



double get_stateSMOLYAK(const int iss, double theta[], const double logcprice, const double qstock, const int ia, const int ieduc){

	double state[DIM_smolyak]; 
		
		int isd = 0; 
		
		state[isd] = exp(theta[0]);  isd++; 
		state[isd] = exp(theta[1]);  isd++; 

		state[isd] = logcprice;    isd++; 
		state[isd] = qstock;    isd++; 
	
	#ifdef DEBUG
	
	if (iss<0 || iss >= DIM_smolyak){ printf("\n Error in get_stateSMOLYAK(): iss = %d \n... exiting...\n", iss); exit(1); }

	
	if (isd != DIM_smolyak) { printf("\n Error in get_stateSMOLYAK(): isd = %d != DIM_smolyak = %d \n ... exiting...\n", isd, DIM_smolyak); exit(1); }

	#endif 

	return state[iss]; 
	
}


int get_stateMODEL(const double state[], double &thetac, double &thetan, double &logcprice, double &qstock, const int ia, const int ieduc){	
	
		int isd = 0; 

		
				thetac = log(state[isd]); isd++;   
				thetan = log(state[isd]); isd++; 
								
				logcprice = (state[isd]); isd++; 								
				qstock    = state[isd]; isd++;

	#ifdef DEBUG			
		if (isd != DIM_smolyak) { printf("\n Error in get_coefEmax(): isd = %d != DIM_smolyak = %d \n ... exiting...\n", isd, DIM_smolyak); exit(1); }
	#endif
	
	return 0; 
	
}


int find_GridIndex(double x, double Mgrid[], int n_mgrid, int &m0, int &m1){
	
		int m=1;  
	
    	if (x >= Mgrid[n_mgrid-1]){   // Above --> linear extrapolation
			
			m0 = n_mgrid-1; 
			m1 = n_mgrid-1-1; 
			
    	}else if (x <= Mgrid[0]){    // Below --> linear extrapolation
			
        	m0 = 0; 
			m1 = 1; 
    	}else {
        	
        	// Find the k's (nearest observed X from left)
        	while(m<n_mgrid && x >= Mgrid[m]){  m++; }
        	
        	m=m-1;  // m is the point in X to the left of Xi
			
			m0 = m; 
			m1 = m+1; 
			
    	} // found m0, m1. 
	
	#ifdef DEBUG
	if ( fabs(Mgrid[m1] - Mgrid[m0]) < 1.0e-10) {  printf("\n Erros in find_GridIndex(), two points too close: mgrid[%d] = %f, and mgrid[%d] = %f \n exitings....\n ", m0, Mgrid[m0], m1, Mgrid[m1]);  }
	#endif 
	
	return 0; 	
		
}



double get_maxV_2choices(const double sige, const double v0, const double v1, const double dv0, const double dv1, double &dmaxV){
// ... this function calculate: max( v0, v1 - sige * eps ), where eps ~ N(0, 1)
	
	double maxV = -1.0e+30; 

	// ... initialization 
	maxV  = v0;
    dmaxV = dv0;
    
                                		
                        if (sige < 1.0e-6){
                            
                            if (v1 > v0){
                                maxV  = v1;
                                dmaxV = dv1;
                            }
                            
                        } else {
                        
						double cut_epsuq = (v1 - v0); 
						double prob    = gsl_cdf_gaussian_P ( cut_epsuq, sige);
						
						if (prob > 1.0e-6){
						
							double e_epsup = cond_mean_cutoff(cut_epsuq, sige); 
							
							maxV  = v0 + prob *  ( v1  - v0 - e_epsup ); 
							dmaxV = dv0 + prob * ( dv1 - dv0 ); 						
//							printf("\n\t\t dmaxV %f, dv0 %f, dv1 %f, prob %f\n", dmaxV, dv0, dv1, prob); 
						}
						
                        } // sige > 1.0e-6						


	return maxV; 
}
